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Executive  Summary 

This  document  represents  the  final  report  of  a  phase  I  STTR  program  between  the  USAF 
Academy  Department  of  Aeronautics  and  Cobalt  Solutions,  LLC.  The  goal  of  the  STTR  is  to  develop  a 
set  of  computational  design  tools  for  closed  loop  flow  control.  These  tools  need  to  cover  the  entire 
design  cycle  from  plant  characterization  through  control  system  development  to  testing  of  the  control 
algorithms  against  truth  models  and  finally  experiments.  The  specific  goal  of  the  phase  I  was  to 
demonstrate  feasibility  of  the  proposed  approach  by  applying  it  to  three  challenging  test  cases:  the 
circular  cylinder  the  “D”  shaped  cylinder,  and  the  NACA  0015  airfoil.  These  three  cases  were  chosen 
because  of  the  differing  approaches  needed  for  each  problem.  Both  cylinders  need  to  be  controlled  by 
stabilizing  the  wake.  The  circular  cylinder  has  variable  separation  points  (depending  on  control  input), 
while  the  “D”  shaped  cylinder  has  fixed  separation.  The  airfoil  may  be  controlled  strictly  be 
controlling  the  separation  (i.e.  eliminating  it).  These  cases  were  all  treated  in  2-D  to  demonstrate 
feasibility  prior  to  attempting  control  in  3-D. 

The  approach  taken  was  to  combine  a  computational  fluid  dynamics  solver  that  is  able  to 
handle  complex  geometries  with  low  dimensional  modeling.  Improvements  were  made  in  the  CFD 
solver  in  order  to  performed  closed  loop  simulations  with  blowing/suction  in  addition  to  the  existing 
motion  capability  and  several  closed  loop  interface  and  post  processing  improvements.  Numerous 
Matlab  routines  were  developed  and  applied  including  a  Proper  Orthogonal  Decomposition  (POD) 
Algorithm,  a  State  Equation  (SEQ)  Algorithm,  a  Linear  Stochastic  Estimation  (LSE)  Algorithm  a 
Closed  Loop  Simulation  Module,  Controllability,  Observability,  and  Stability  (COS)  Algorithms,  and 
Sensor  Placement  and  Mode  Estimation  Modules  (SME).  The  toolbox  was  applied  to  the  three  cases 
outlined  with  demonstrated  success  on  all  three.  The  cylinder  wakes  were  stabilized  to  a  large  extent 
with  a  reduction  in  lift  oscillations  and  drag.  The  separation  on  the  NACA  0015  was  reduced,  though 
not  eliminated  (in  the  mean).  The  tools  were  very  useful  in  designing  controllers  for  all  three  cases. 
Conclusions  and  lessons  learned  from  these  test  cases  are  outlined  with  recommendations  made  for 
proceeding  to  phase  II  of  the  STTR. 
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1.  Introduction 

1.1  Significance  of  Problem 

The  main  idea  of  flow  control  is  the  improvement  of  aerodynamic  characteristics  of  air  vehicles 
and  munitions  enabling  augmented  mission  performance.  Flow  control  can  either  be  passive  or  active 
depending  on  whether  energy  is  added  to  the  flow.  Furthermore,  active  flow  control  may  be 
characterized  by  open-loop  or  closed-loop  techniques.  Gad-el-Hak  (1996)  provides  an  insight  into  the 
advances  in  the  field  of  flow  control.  Research  into  closed-loop  flow  control  methods  has  increased 
over  the  past  two  decades.  Cattafesta  et  al  (2003)  provide  a  useful  classification  of  active  flow  control 
and  describes  the  main  components  of  a  feedback  control  system.  Before  proceeding  into  the  details  of 
modeling  and  control,  it  is  imperative  to  appreciate  the  reasons  as  to  why  closed-loop  control  is  of 
consequence  anti  the  main  advantages  associated  with  its  implementation  in  flow  control  problems.  It 
is  advantageous  to  opt  for  closed-loop  flow  control  for  it: 

•  Enables  addressing  problems  that  have  over  the  years  not  been  solved  using  passive  means  and  /or 
open-loop  techniques. 

•  Provides  performance  augmentation  of  an  open-loop  flow  control  system. 

•  Lowers  the  amount  of  energy  required  to  manipulate  the  flow  to  induce  the  desired  behavior.  This 
aspect  affects  actuation  requirement  and  may  be  a  deciding  factor. 

•  Enables  adaptability  to  a  wider  operating  envelope,  thereby  limiting  the  drop  in  performance 
associated  with  multiple  design  working  points. 

•  Provides  design  flexibility  and  robustness. 

Feedback  flow  control  is  an  emerging  discipline  which  aims  at  applying  the  methods  of  closed 
loop  feedback  control  to  problems  in  fluid  dynamics.  Traditionally,  these  two  research  fields  have 
existed  independently  and  without  mutual  involvement.  Control  scientists  have  applied  their 
knowledge  historically  to  problems  with  a  relatively  small  dimensionality.  Single  input  single  output 
systems  of  linear  behavior  are  very  well  understood  and  theoretically  covered.  More  recently,  feedback 
control  has  been  applied  to  more  complex  problems  in  structural  dynamics  as  well  as  fluid-structure 
interaction.  These  problems  are  of  higher  dimensionality  and  typically  can  be  described  by  a  small 
number  of  dominant  modes.  They  also  often  involve  multiple  sensors  and  actuators.  Feedback  control 
of  relatively  complex  structures  is  currently  state  of  the  art  in  control  sciences,  and  very  few  control 
researchers  have  ventured  into  problems  of  higher  complexity.  Several  applications  of  closed-loop 
control  have  been  reported  in  literature.  Specific  areas  of  interest  include  flow-induced  cavity 
resonance  (Cattafesta,  2003),  vectoring  control  of  a  turbulent  jet,  separation  control  of  the  NASA 
Langley  hump  model  (a  variation  on  the  Glauert  Glas  II  airfoil)  control  of  the  vortex  motion  in  the 
combustion  recirculation  region  and  control  of  vortex  shedding  in  circular  cylinder  wakes. 

Even  simple  problems  in  fluid  dynamics,  however,  present  themselves  to  the  control  scientist 
as  a  system  of  much  higher  dimensionality  whose  governing  equations  are  highly  nonlinear  and,  for 
most  all  boundary  conditions  of  technical  interest,  cannot  be  solved  in  closed  form.  This  poses  a 
formidable  challenge  for  the  application  of  feedback  flow  control  to  fluid  dynamics  problems  of 
technical  interest.  On  the  other  hand,  since  most  fluid  dynamics  problems  are  featuring  local  or  global 
instabilities,  feedback  flow  control  holds  great  promise  in  improving  these  flow  fields  with  very 
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moderate  control  inputs,  while  yielding  large  reductions  in  flow  quantities  of  technical  importance  like 
drag  and  unsteadiness  of  loads  imposed  by  the  flow  field  onto  structures  exposed  to  the  flow. 

The  biggest  challenge  to  be  overcome  in  order  to  be  able  to  make  progress  in  the  field  of 
feedback  flow  control  is  in  the  area  of  modeling  the  dynamics  of  the  flow  field.  So  far,  there  have  been 
three  different  approaches  to  modeling.  Fluid  dynamics  researchers  usually  refrained  from  modeling 
the  flow  at  all.  Instead,  they  relied  on  empirical  tuning  of  black  box  type  controllers.  Their  attempt  to 
control  flow  fields  was  met  with  very  moderate  success,  and  only  for  very  simple  flow  fields  were  they 
able  to  achieve  significant  improvements.  The  linear  control  strategies  applied  were  typically  limited  to 
single  sensor  /  single  actuator  scenarios  and  important  questions  of  controllability  and  observability 
could  not  be  addressed  due  to  the  nonexistence  of  a  mathematical  model  which  is  a  prerequisite  for  this 
kind  of  investigations. 

The  second  approach,  usually  employed  by  researchers  with  a  background  in  applied 
mathematics,  is  based  on  solving  the  Navier  Stokes  equations  and  their  respective  boundary  conditions 
using  discretization  on  a  grid  and  numerical  methods.  While  this  approach  is  shared  by  the  field  of 
computational  fluid  dynamics  (CFD),  for  application  to  feedback  control  not  just  the  forward  problem 
that  CFD  solves  but  also  the  inverse  problem  needs  to  be  considered  especially  when  optimum  controls 
approaches  are  considered.  This  leads  to  various  forms  of  the  Ricatti  equations,  which  have  proven  to 
be  much  more  difficult  to  solve  than  the  Navier  Stokes  equations  themselves.  For  most  fluid  dynamics 
problems  a  solution  of  the  Ricatti  equations  with  current  state  of  the  art  computing  methods  seems  out 
of  reach.  The  problems  encountered  by  both  of  these  approaches  to  feedback  flow  control  led  to  the 
third  approach  which  employs  low  dimensional  modeling.  This  technique  has  been  applied  with  great 
success  to  control  problems  involving  structural  dynamics.  The  goal  is  to  capture  the  relevant 
dynamics  of  the  system,  be  it  a  mechanical  structure  or  a  flow  field,  by  a  finite  number  of  spatial 
modes.  The  state  of  the  system  is  then  characterized  by  these  modes,  and  an  observer  or  state  estimator 
is  usually  employed  to  derive  the  state  of  the  system  in  terms  of  the  mode  amplitudes  based  on  sensor 
readings.  The  controller  then  acts  on  the  mode  estimates  as  opposed  to  the  sensor  readings.  The 
advantage  of  this  approach  is  in  the  ability  to  actually  model  the  important  features  of  the  system 
without  having  to  solve  the  governing  equations  directly  or  numerically.  Furthermore,  the  model  can 
be  derived  from  either  experimental  or  simulation  data,  thus  allowing  development  of  a  controller 
based  on  the  dynamics  of  the  actual  system  that  is  to  be  controlled.  Of  the  three  described  attempts  at 
feedback  flow  control,  the  use  of  low  dimensional  modeling  has  shown  the  greatest  promise  to  date. 

1.2  Main  Objectives 

This  STTR  research  program  aims  at  developing  computational  tools  necessary  to  develop 
feedback  control  algorithms  based  on  low  dimensional  models.  These  tools  need  to,  cover  the  entire 
design  cycle  from  plant  characterization  through  control  system  development  to  testing  of  the  control 
algorithms  against  truth  models  and  finally  experiments.  Currently,  no  such  tools  are  commercially 
available,  and  researchers  interested  in  feedback  flow  control  have  to  develop  their  own  tools  from 
scratch.  Since  neither  controls  nor  fluids  experts  are  typically  expert  programmers,  they  need  to  expend 
significant  time  and  effort  to  even  obtain  tools  of  limited  functionality,  typically  only  usable  for  a  very 
specific  problem.  The  main  goal  of  this  STTR  collaboration  between  the  US  Air  Force  Academy  and 
Cobalt  Solutions  is  to  join  the  controls  and  fluids  expertise  at  the  Academy  with  the  programming  and 
numerical  mathematical  solution  skills  at  Cobalt  Solutions.  The  result  will  be  a  toolbox  that  provides 
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an  integrated  closed  loop  flow  control  development  environment  covering  the  entire  range  of  tasks 
from  fluid  dynamics  simulation  through  low  order  modeling,  mathematical  model  development  and 
controller  design  to  controller  implementation  and  testing  of  the  controller  design  by  applying  it  to  a 
truth  model  based  on  closed  loop  fluid  dynamic  simulation. 

The  toolbox  enables  closed  loop  CFD  simulations  based  on  Low-Dimensional  Models  obtained 
from  POD.  The  toolbox  features  the  following  main  components: 

■  Computational  Fluid  Dynamics  (CFD)  Solver 

■  Proper  Orthogonal  Decomposition  (POD)  Algorithm 

■  State  Equation  (SEQ)  Algorithm 

■  Linear  Stochastic  Estimation  (LSE)  Algorithm 

■  Closed  Loop  Simulation  Module 

■  Controllability,  Observability,  and  Stability  (COS)  Algorithms 

•  Sensor  Placement  and  Mode  Estimation  Modules  (SME) 

These  tools  cover  the  entire  controller  design  cycle  from  model  development  through  control 
algorithm  design  to  verification  of  the  design  on  the  truth  model. 

1.3  Toolbox  Application 

In  order  to  demonstrate  and  evaluate  the  performance  of  the  toolbox  components  developed,  we 
plan  to  use  three  (3)  different  benchmark  problems: 

•  Circular  Cylinder 

•  D-  Shaped  Cylinder 

•  Symmetric  airfoil  (NACA  00 1 5)  at  high  angle  of  attack 

While  the  first  two  geometries  develop  a  typical  wake  flow  featuring  an  absolute  instability,  the 
third  geometry  aims  at  separation  control. 

Absolutely  unstable  wake  flows  have  proven  to  be  one  of  the  most  intractable  problems  in  fluid 
dynamics  using  active  open  loop  flow  control.  Typically,  the  unsteadiness  of  the  flow  is  actually 
increased  along  with  an  increase  in  drag,  or  the  small  performance  gain  achieved  using  high  amplitude 
forcing  at  high  frequencies  does  not  outweigh  the  energy  expended  for  the  forcing.  Thus  this  type  of 
flow  is  an  ideal  benchmark  problem  for  feedback  flow  control,  and  success  at  very  small  Reynolds 
numbers  has  been  achieved. 

While  separation  control  has  been  successfully  addressed  in  literature  using  open  loop  forcing 
techniques,  it  will  be  interesting  to  see  what  performance  gain  can  be  achieved  using  closed  loop 
techniques.  This  makes  the  airfoil  at  high  angle  of  attack  an  interesting  feedback  flow  control  problem. 

The  different  forcing  methods  outlined  in  the  previous  section  will  be  applied  to  demonstrate 
the  capabilities  of  the  system.  State  estimation  will  be  obtained  both  using  sensors  distributed  in  the 
flow  field,  as  well  as  surface  mounted  sensors.  Linear  stochastic  estimation  will  be  employed  for  the 
surface  mounted  sensors  to  map  their  readings  to  POD  mode  amplitude  estimates.  For  all  three 
geometries,  we  will  keep  the  Reynolds  numbers  low  enough  to  be  able  to  obtain  physical  results 
without  the  use  of  turbulence  models. 
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This  report  paper  is  organized  as  follows:  The  next  section  provides  an  overview  on  the 
computational  tools  developed  within  the  framework  of  this  Phase  I  STTR  program.  The  modeling  and 
control  aspects  are  presented  in  the  following  section,  and  the  applications  of  the  design  tools  are 
presented  in  the  following  three  chapters.  The  first  application  is  the  circular  cylinder  followed  by  the 
D-  shaped  cylinder  and  finally  by  a  symmetric  airfoil  (NACA  0015).  Then,  a  chapter  summarizes 
Phase  I  and  presents  recommendations  for  Phase  II.  At  the  end  of  the  report,  we  have  included  a  listing 
of  the  main  Matlab  files  developed  in  the  areas  of  modeling,  estimation  and  control. 
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2.  Computational  Tools  Development 

2.1  Overview  of  COBALT  CFD  Solver 

The  CFD  solver  used  is  Cobalt,  a  commercial  unstructured  finite-volume  code  developed  for 
solution  of  the  compressible  Navier-Stokes  equations.  The  basic  algorithm  is  described  in  Strang  et  al. 
(1999),  although  substantial  changes/improvements  have  been  made  since  this  paper.  The  numerical 
method  is  a  cell-centered  finite  volume  approach  applicable  to  arbitrary  cell  topologies  (e.g, 
hexahedrals,  prisms,  tetrahdra).  The  spatial  operator  uses  the  exact  Riemann  Solver  of  Gottlieb  and 
Groth  (1988),  least  squares  gradient  calculations  using  QR  factorization  to  provide  second  order 
accuracy  in  space,  and  TVD  flux  limiters  to  limit  extremes  at  cell  faces.  A  point  implicit  method  using 
analytic  first-ortjer  inviscid  and  viscous  Jacobians  is  used  for  advancement  of  the  discretized  system. 
For  time-accurate  computations,  a  Newton  sub-iteration  scheme  is  employed,  the  method  is  second 
order  accurate  in  time. 

For  parallel  performance,  Cobalt  utilizes  the  domain  decomposition  library  ParMETIS 
(Karypis  et  al.  1997)  to  provide  optimal  load  balancing  with  a  minimal  surface  interface  between 
zones.  Communication  between  processors  is  achieved  using  Message  Passing  Interface  (MPI),  with 
parallel  efficiencies  above  95%  on  as  many  as  1024  processors  (Grismer  et  al.  1998). 

The  vast  majority  of  Air  Force  vehicles  operate  at  high  Reynolds  numbers  where  the  fluid  is 
turbulent.  The  main  methods  for  calculating  turbulent  flows  with  a  CFD  solver  are  Direct  Numerical 
Simulation  (DNS),  Large  Eddy  Simulation  (LES),  and  Reynolds-averaged  Navier  Stokes  (RANS). 
The  RANS  approach  is  the  most  economical  since  it  is  designed  to  solve  for  the  mean  flow,  but  is 
subject  to  many  modeling  approximations.  Since  it  models  rather  than  resolves  the  bulk  if  not  all  of 
the  turbulent  motions  this  is  an  inappropriate  choice  for  most  flow  control  applications.  DNS,  on  the 
other  hand,  makes  no  modeling  assumption  but  is  the  most  expensive  approach  since  all  turbulent 
motions  must  be  resolved  by  the  grid.  Since  the  smallest  scales  of  turbulence  (the  Kolmogorov  length 
scale)  decrease  with  Reynolds  number,  this  limits  this  approach  to  low  Reynolds  number  flows.  LES 
is  less  expensive  than  DNS  since  it  models  only  the  small  subgrid  scales  of  motion  and  resolves  the 
rest  of  the  turbulent  motions.  However,  since  the  “large”  scales  in  the  boundary  layer  are  on  the  order 
of  the  boundary  layer  thickness  (which  is  quite  thin  for  high  Reynolds  number  flows),  this  method  is 
cost  prohibitive  at  high  Reynolds  numbers  for  wall  bounded  flows.  ' 

Detached-Eddy  Simulation  (DES)  is  a  hybrid  technique  first  proposed  by  Spalart  et  al  (1997) 
for  prediction  of  turbulent  flows  at  high  Reynolds  numbers  (see  also  Spalart  2000).  The  motivation  for 
developing  DES  was  to  combine  the  most  favorable  aspects  of  RANS  and  LES,  i.e.,  the  acceptable 
predictions  using  RANS  models  of  thin  shear  layers  (e.g.,  attached  boundary  layers)  and  LES  for 
resolution  of  time-dependent,  three-dimensional  large  eddies.  For  natural  applications  of  DES,  RANS 
is  applied  in  the  boundary  layer,  while  outside  the  boundary  layer  in  the  separated  region.  An  array  of 
flows  ranging  from  “building  block”  applications  such  as  the  flow  over  a  cylinder,  sphere,  aircraft 
forebody,  and  missile  base  to  complex  geometries  including  full,  aircraft  have  been  modeled  (e.g.,  see 
Travin  et  al.  2001,  Squires  et  al.  2001,  Constantinescu  et  al.  2002,  Forsythe  et  al.  2002,  Hansen  and 
Forsythe  2003,  Constantinescu  et  al.  2003).  These  and  other  applications  have  been  largely  successful, 
illustrating  DES  capabilities  in  accurately  resolving  chaotic  unsteady  features  in  the  separated  regions 
along  with  a  rational  treatment  of  the  attached  boundary  layers.  Recent  DES  predictions  of  the  flow 
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around  complex  configurations  (all  with  Cobalt)  include  the  massively  separated  flow  around  an  F- 
15E  at  65°  angle  of  attack  reported  by  Forsythe  et  al.  (2004)  (the  first  eddy-resolving  simulation  of 
flow  around  a  full  aircraft  configuration),  transonic  shock-separated  flow  over  an  F/A-18E  by  Forsythe 
and  Woodson  (2003),  and  vortex  breakdown  on  an  F-l  8C  by  Morton  et  al.  (2003). 

The  use  of  a  pre-existing  CFD  code  that  has  been  thoroughly  verified  within  the  context 
of  turbulence  resolving  calculations  has  reduce  the  costs/risks  of  the  proposed  work.  The  Phase  I  work 
with  the  CFD  solver  has  consisted  of  integrating  the  solver  more  tightly  with  the  flow  control  tools  and 
building  a  3D  baseline  of  open  loop  applications. 

2.2  Development  and  Coding  Issues 

First,  it  was  decided  that  a  Hierarchical  Data  Format  (HDF)  file  would  best  serve  as  the 
communication?  interface  between  the  fluid  flow  solver,  Cobalt,  and  the  USAFA-developed  feedback 
controller.  HDF,  developed  by  the  National  Center  for  Supercomputing  Applications,  is  a  multi-object 
file  format  for  sharing  scientific  data  in  a  distributed  environment.  HDF  is  a  standard  file  format  in  the 
science  and  engineering  field,  facilitating  the  incorporation  of  feedback  controllers  developed  outside 
of  USAFA.  Also,  the  platform  independence  of  HDF  provides  considerable  flexibility  in  operating  the 
feedback-controller-Cobalt-HDF  system.  For  example.  Cobalt  could  be  running  on  128-processor 
parallel  computer,  with  the  feedback  controller  running  on  an  individual’s  workstation  and  the  HDF 
file  residing  on  a  fileserver.  The  sole  requirement  is  that  the  three  machines  are  networked  together. 

Second,  it  was  agreed  upon  that  the  feedback  controller  would  control  two  similar  boundary 
conditions  in  the  fluid  flow  solver.  Both  boundary  conditions  simulate  flow  entering  or  leaving  the 
fluid  domain;  differences  between  the  two  lie  in  the  details.  In  the  first  method,  called  the  “mass  flow 
boundary  condition”,  the  external  controller  controls  the  mass  flow  entering  or  leaving  the  domain  as 
well  as  the  total  temperature  and  total  pressure  of  that  mass  flow.  This  boundary  condition  is  valid  for 
incompressible  and  compressible  flows.  The  second  boundary  condition,  called  the  “velocity  boundary 
condition”,  is  based  on  a  method  developed  by  S.  A.  Morton  of  USAFA.  This  method  is  strictly  valid 
for  incompressible  flows  only,  but  can  be  used  with  good  results  up  to  a  Mach  number  of  around  0.2. 
Here,  the  external  controller  controls  only  the  velocity  of  the  fluid  entering  or  leaving  the  domain; 
pressure  and  temperature  conditions  can  be  derived  based  on  the  knowledge  of  the  temporal  derivative 
of  velocity  (calculated  either  by  the  external  controller  or  by  the  fluid  flow  solver)  with  the  assumption 
of  incompressible  flow.  , 

The  mass  flow  boundary  condition  was  first  added  to  Cobalt  and  validated.  The  validation  was 
straightforward,  since  this  boundary  condition  is  the  simple  union  of  two  existing,  and  well  validated, 
boundary  conditions  in  Cobalt.  The  coding  required  for  Cobalt  to  create,  read  from,  and  write  to  an 
HDF  file  was  then  accomplished.  A  very  simple  external  controller  was  then  developed  to  validate  that 
information  was  being  properly  passed  via  the  HDF  file  for  a  simple  test  case.  Results  of  this  test  were 
compared  with  those  resulting  from  a  similar  test  where  Cobalt  itself  varied  the  mass  flow  boundary 
condition  in  exactly  the  same  manner  as  the  external  controller.  Results  between  the  two  tests  were 
identical  to  with  machine  round-off  error. 

Coding  of  the  velocity  boundary  condition  was  accomplished.  Since  this  boundary  condition  is 
new  to  Cobalt,  it  was  subjected  to  extensive  validation.  However,  once  this  validation  was  complete, 
integrating  it  with  the  HDF  file  interface  advanced  quickly,  since  much  of  the  integration  for  the  mass 
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flow  boundary  condition  was  reused.  During  the  coding,  an  “open  -loop”  boundary  condition  was  also 
added  that  did  not  require  interfacing  through  HDF  to  an  external  controller.  This  module  was  used  for 
the  prediction  of  the  flow  over  a  wall  mounted  hump  by  Krishnan  et  al.  (2004). 


Fig.  2.1:  Open  loop  blowing/suction  using  on  the  wall  mounted  hump  by  Krishnan  et  al.  using  Cobalt. 

Contours  of  streamwise  velocity  with  instantaneous  streamlines. 

Cobalt  Solutions  provided  data  exchange  interfaces  in  the  CFD  solver  that  allowed  both 
unforced  and  open  loop  forced  data  collection  to  HDF  format  files,  as  well  as  control  interaction 
during  a  CFD  simulation.  The  latter  feature  allowed  the  CFD  Solver  (Cobalt)  to  be  used  as  a  truth 
model  to  verify  a  given  controller  design.  By  using  files  in  order  to  interface  the  CFD  solver  and  the 
control  toolbox,  debugging  of  the  overall  system  was  facilitated.  Additionally,  the  entire  CFD  code 
was  not  needed  to  be  recompiled  when  new  control  algorithms  were  implemented.  Within  the  Phase  I, 
the  following  limitations  were  applied: 

•  All  modules  other  than  the  CFD  Solver  were  implemented  in  a  high  level  language,  in 
particular,  in  Matlab.  Porting  of  Matlab  algorithms  to  lower  level  languages  and  parallelization 
of  these  codes  will  be  performed  in  Phase  II,  in  order  to  accommodate  larger  (3D)  data  sets  and 
improve  performance. 

•  All  closed  loop  simulations  were  two-dimensional  in  order  to  keep  the  required  computation 
times  low,  while  still  being  able  to  achieve  the  objective  of  demonstrating  the 'performance  of 
the  tools.  While  open  loop  three-dimensional  CFD  performance  was  demonstrated  in  Phase  I, 
closed  loop  3D  simulations  will  be  performed  in  Phase  II. 

•  All  of  the  Matlab  tools  were  implemented  for  two-dimensional  datasets.  Extension  of  the  tool’s 

capabilities  to  three-dimensions  will  take  place  in  Phase  II.  « 

In  order  to  perform  closed  loop  flow  control,  actuation  is  needed.  The  following  actuation 
approaches  was  developed  and  validated: 

•  Translation  normal  to  the  flow 

•  Blowing  and  suction  at  the  model  surface  ' 
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The  toolbox  allows  for  an  unlimited  number  of  sensors  to  be  placed  within  the  computational  domain. 
At  any  sensing  location,  the  following  flow  variables  will  be  available  to  the  controller: 

•  All  three  velocity  components  (U,  V,  W) 

•  Pressure 

•  Density 

Additionally,  the  following  global  flow  properties  are  available  to  the  controller  and  can  be  used  for 
closing  the  feedback  loop: 

•  X,Y,Z  components  of  force  acting  on  the  body 

•  X,Y,Z  components  of  moment 

As  a  result  of  the  current  research  program,  Cobalt  can  now  output  instantaneous  flow  data  at 
arbitrary  user-sj^cified  locations,  ‘taps’,  in  Hierarchical  Data  Format  (HDF).  Previously,  only  ASCII 
format  was  available.  Due  to  the  nature  of  HDF,  all  flow  quantities  available  for  post-processing,  16  in 
all,  are  output  for  every  tap.  Optimization  of  the  algorithm  resulted  in  a  very  low  interpolation 
overhead  during  the  output  of  only  1 5%  of  the  computational  time  per  time  step,  for  a  large  number  of 
taps.  Note  that  this  overhead  is  only  incurred  during  data  output.  Furthermore,  the  implementation  of 
HDF  output  for  an  external  controller,  for  both  rigid  body  motion  and  blowing  and  suction  boundary 
patches  was  successfully  completed.  An  additional  improvement  is  that  the  flow  quantities  are  now  tri- 
linearly  interpolated  to  the  taps.  Previously,  the  computational  cell  containing  any  given  tap  was  first 
found  and  the  flow  state  at  that  cell  centroid  was  then  assigned  to  the  tap  location.  The  current  method, 
while  more  accurate,  requires  additional  CPU  time  over  the  original  method.  These  two  modifications 
have  enabled  more  rapid  and  more  accurate  POD  creation. 

2.3  Debugging  Issues 

After  development  and  coding  of  the  toolbox  components  was  accomplished,  a  thorough 
debugging  phase  followed.  We  compared  the  results  of  the  toolbox  modules  to  published  literature 
wherever  possible  and  available.  Test  cases  were  run  to  first  verify  the  performance  of  the  individual 
toolbox  components  separately.  After  the  individual  toolbox  components  performance  were  verified, 
the  integration  of  the  toolbox  components  with  each  other  was  tested,  with  particular  attention  on  the 
interface  between  the  CFD  solver  and  the  Matlab  based  components  of  the  toolbox. 

Finally,  the  external  controller  HDF  interface  was  debugged  in  collaboration  with  ySAFA. 
After  debugging  was  completed,  we  applied  the  toolbox  to  the  benchmark  problems  outlined  in 
Section  1.3  of  this  report. 

2.4  Computational  Highlights 

Overall,  the  final  toolbox  will  be  able  to  handle  any  flow  geometry  at  any  Reynolds  number. 
The  CFD  solver  is  capable  of  handling  unstructured  grids.  It  can  perform  simulations  from  Euler 
through  laminar  Navier  Stokes  up  to  a  variety  of  turbulence  models  as  well  as  Detached-Eddy 
Simulation  (DES).  Rigid  body  motion  as  well  as  variable  boundary  conditions  are  implemented.  On 
the  toolbox  side,  all  of  the  above  features  are  supported,  while  limiting  the  flow  geometry  to  two- 
dimensional  at  this  point  in  time. 
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One  of  the  main  integration  goals  reached  was  to  interface  the  toolbox  to  the  CFD  simulation 
using  open  standard  data  formats,  in  particular,  the  Hierarchical  Data  Format  (HDF)  developed  by  the 
National  Center  for  Supercomputing  Applications  (NCSA).  The  open  standard  data  format  will  allow 
users  to  use  their  plotting  and  data  post  processing  software  tools  of  choice,  since  most  major  software 
packages  support  this  format.  Additionally,  this  data  format  is  very  efficient  in  terms  of  file  size  while 
providing  convenient  and  quick  access  to  the  data. 
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3.  Modeling,  Estimation  and  Control  Tools  Development 

3.1  Modeling  Approaches 

A  closed-loop  flow  control  system  is  comprised  of  a  controller  that  introduces  a  perturbation 
into  the  flow,  via  a  set  of  actuators,  to  obtain  desired  performance.  Furthermore,  the  controller  acts 
upon  information  provided  by  a  set  of  sensors.  There  are  three  basic  approaches  to  closed-loop  flow 
control  of  a  two-dimensional  wake  as  depicted  in  Figure  3.1 . 

Model  Independent  Approach  -  involves  the  introduction  of  sensors  in  the  wake  and  using  a  control 
law  (usually  linear)  which  produces  a  command  to  the  actuator  that  forces  the  flow.  The  advantage  of 

this  approach  is^o  show: 

* 

•  No  model  of  the  flow  field  is  required  for  controller  design. 

•  Direct  feedback  eliminates  the  need  for  a  state  estimator. 

•  A  simple  control  law  may  be  implemented  in  an  experimental  set  up  with  relative  ease. 

For  example,  let  us  consider  the  circular  cylinder  wake  problem.  Experimental  studies  show 
that  a  linear  proportional  feedback  control  based  on  a  single  sensor  feedback  is  able  to  delay  the  onset 
of  the  wake  instability,  rendering  the  wake  stable  at  Re  about  20%  higher  than  the  unforced  case. 
Above  Re  =  60,  a  single-sensor  feedback  may  suppress  the  original  mode  but  destabilizes  one  of  the 
other  modes  (Roussopoulos,  1993).  This  approach  is  relatively  simple  to  implement  experimentally. 
However,  the  results  are  very  limited  for  the  challenging  problem  of  an  absolutely  unstable  wake. 


•  Simple  to  •  Ideal  control  •  Recent  developments 


implement 

experimentally 


approach,  complete  in  effective  low- 

set  of  equations  dimensional  models 


•  Works  for  a 
limited  number 
of  applications 


•  Computationally 
intensive 

•  Can  be  implemented 
in  restricted  cases 


•  Can  be  implemented 
with  relative  ease 

•Model  building  is 
tough 


Fig.  3.1:  Approaches  to  Closed-Loop  Flow  Control 
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Direct  Navier  Stokes  Approach  -  This  approach  is  more  structured  as  it  applies  conventional  and 
proven  model-based  control  strategies  such  as  optimal  control  theory  for  flow  control  problems. 
Abergel  and  Temam  (1990)  developed  conditions  for  optimality  for  a  few  simple  applications. 
However,  real  time  implementation  of  this  approach  to  the  cumbersome  unsteady  Navier-Stokes 
equations  is  not  practical. 

Low-Dimensional  Approach  -  Low-dimensional  modeling  is  a  vital  building  block  when  it  comes  to 
realizing  a  structured  model-based  closed-loop  strategy  for  flow  control.  For  control  purposes,  a 
practical  procedure  is  needed  to  break  down  the  velocity  field,  governed  by  Navier  Stokes  partial 
differential  equations,  by  separating  space  and  time.  A  common  method  used  to  substantially  reduce 
the  order  of  the^nodel  is  proper  orthogonal  decomposition  (POD).  This  method  is  an  optimal  approach 
in  that  it  will  capture  a  larger  amount  of  the  flow  energy  in  the  fewest  modes  of  any  decomposition  of 
the  flow.  The  POD  method  may  be  used  to  identify  the  characteristic  features,  or  modes,  of  a  cylinder 
wake  as  demonstrated  by  Gillies  (1998). 

The  major  building  blocks  of  this  structured  approach  are  comprised  of  a  reduced-order  POD 
model,  a  state  estimator  and  a  controller.  The  desired  POD  model  contains  an  adequate  number  of 
modes  to  enable  reasonable  modeling  of  the  temporal  and  spatial  characteristics  of  the  large  scale 
coherent  structures  inherent  in  the  flow  though  it  may  faithfully  reproduce  the  flow.  Further  details  of 
the  POD  method  may  be  found  in  the  book  by  Holmes,  Lumley,  and  Berkooz  (1996).  A  common 
approach  referred  to  as  the  method  of  “snapshots”  introduced  by  Sirovich(1987)  is  employed  to 
generate  the  basis  functions  of  the  POD  spatial  modes  from  flow-field  information  obtained  using 
either  experiments  or  numerical  simulations.  This  approach  to  controlling  the  global  wake  behavior 
behind  a  circular  cylinder  was  effectively  employed  by  Gillies  (1998)  and  is  also  the  approach 
followed  in  this  research  effort. 

For  low-dimensional  control  schemes  to  be  implemented,  a  real-time  estimation  of  the  modes 
present  in  the  wake  is  necessary  since  it  is  not  possible  to  measure  them  directly  especially  in  real¬ 
time.  An  illustration  of  the  various  blocks  within  the  low-dimensional  modeling  approach  is  presented 
in  Fig.  3.2.  Velocity  field  data,  provided  from  either  simulation  or  experiment,  is  fed  into  the  POD 
procedure.  The  time  histories  of  the  temporal  coefficients  of  the  POD  model  are  determined  by 
applying  the  least  squares  technique  to  the  spatial  Eigen-functions  and  the  unforced  flow.  Then,  the 
estimation  of  the  low-dimensional  states  is  provided  using  a  linear  stochastic  estimator  (LSE).  Sensor 
measurements  may  take  the  form  of  wake  velocity  measurements,  as  in  this  effort,  or  body-mounted 
pressure  measurements.  This  process  leads  to  the  state  and  measurement  equations,  required  for  design 
of  the  control  system.  For  practical  applications  it  is  desirable  to  reduce  the  information  required  for 
estimation  to  the  minimum.  The  requirement  for  the  estimation  scheme  then  is  to  behave  as  a  modal 
filter  that  “combs  out”  the  higher  modes.  The  main  aim  of  this  approach  is  to  thereby  circumvent  the 
destabilizing  effects  of  observation  “spillover”  as  described  by  Balas  (1978).  Spillover  has  been  the 
cause  for  instability  in  the  control  of  flexible  structures  and  modal  filtering  was  found  to  be  an 
effective  remedy. 
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Basic  Modeling  Approach 


Fig.  3.2:  Low-Dimensional  Modeling  Strategy 


The  intention  of  the  proposed  strategy  is  that  the  signals,  provided  by  a  certain  configuration  of 
sensors  placed  in  the  wake,  are  processed  by  the  estimator  to  provide  the  estimates  of  the  first  two 
modes.  The  estimation  scheme,  based  on  the  linear  stochastic  estimation  procedure  introduced  by 
Adrian  (1997),  predicts  the  temporal  amplitudes  of  the  first  two  POD  modes  from  a  finite  set  of 
measurements  obtained  from  either  computational  or  experimental  data.  A  major  design  challenge  lies 
in  finding  an  appropriate  number  of  sensors  and  their  locations  that  will  best  enable  the  desired  modal 
filtering. 

3.2  Proper  Orthogonal  Decomposition  (POD) 

Feasible  real  time  estimation  and  control  of  the  cylinder  wake  may  be  effectively  realized  by 
reducing  the  model  complexity  of  the  cylinder  wake  as  described  by  the  Navier-Stokes  equations, 
using  POD  techniques.  POD,  a  non-linear  model  reduction  approach  is  referred  to  in  the  literature  as 
the  Karhunen-Loeve  expansion.  The  desired  POD  model  will  contain  an  adequate  number  of  modes  to 
enable  modeling  of  the  temporal  and  spatial  characteristics  of  the  large-scale  coherent  structures 
inherent  in  the  flow,  but  no  more  modes  than  necessary. 
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The  basis  functions  of  the  POD  spatial  modes  may  be  obtained  from  the  numerical 
solution  of  the  Navier-Stokes  equations  using  CFD  simulations  or  by  incorporating  experimental  data. 
For  control  design  purposes,  the  POD  method  enables  the  Navier-Stokes  equations  to  be  modeled  as  a 
set  of  ordinary  differential  equations  (O.D.E.).  The  decomposition  of  this  component  of  the  velocity 
field  is  as  follows: 


u(x,y,t)  =  U(x,y)  +  u(x,y,t )  (3.1) 

where  (/[m/s]  denotes  the  mean  flow  and  w[m/s]  is  the  fluctuating  component  that  may  be  expanded  as: 

n 

u(x,  y,  0  =  2  ak  (t)0-k}  (x,  y)  (3.2) 

k= 1 

where  o*(t)  denotes  the  time-dependent  coefficients  and  (p,(x,y)  represents  the  non-dimensional  spatial 
Eigenfunctions  determined  from  the  POD  procedure.  In  some  cases  the  velocity  component  may 
represent  the  vorticity,  which  is  calculated  from  the  flow  field. 

Next,  the  empirical  correlation  matrix,  R  is  computed.  A  simple  approximation  technique  is 
applied  to  obtain  the  numerical  integration.  In  this  effort,  the  correlation  matrix  is  computed  using  the 
inner  product.  Solving  the  Eigenvalue  problem,  the  Eigenvalues  and  the  orthogonal  spatial 
Eigenfunctions,  <pi(x,y)  are  obtained.  Since  the  Eigenvalues  measure  the  relative  energy  of  the  system 
dynamics  contained  in  that  particular  mode,  they  may  be  normalized  to  correspond  to  a  percentage. 
Finally,  the  time  histories  of  the  temporal  coefficients  of  the  POD  model,  ak(t),  are  determined  using 
the  extracted  spatial  modes  and  the  data  of  the  unforced  flow.  For  an  arbitrarily  forced  circular 
cylinder,  we  can  write  the  low-dimensional  wake  model  as: 

da.  , 

~^~  =  Sk(an)  +  bkfa  (3.3) 

t 

where  gk,  for  k  modes,  is  a  quadratic  non-linear  function  of  the  time-dependent  mode  coefficients.  bk 
are  coefficients  associated  with  the  control  input  and  fa  is  the  feedback  control  input  to  the  cylinder. 
For  the  open-loop  case  fa  =  0.  For  a  full  state  feedback  system,  the  closed  loop  control  input,  fa,  is  a 
function  of  ak.  However,  it  is  not  possible  to  obtain  a  direct  measurement  of  ak.  An  essential  aspect  of 
reduced  order  modeling  concerns  truncation.  How  many  modes  are  required  and  what  are  the  criteria 
for  effective  truncation?  The  answers  to  this  question  has  been  partially  addressed  by  Cohen  et  al 
(2003)  for  the  example  of  the  circular  cylinder  wake  control.  This  effort  showed  that  control  of  the 
POD  model  of  the  von  Karman  vortex  street  in  the  wake  of  a  circular  cylinder  at  Re  =  100  is 
enabled  using  just  the  first  mode.  - 
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Furthermore,  feedback  based  on  the  first  mode  alone  suppressed  all  the  other  modes  in 
the  four  mode  POD  model.  At  this  point,  it  is  imperative  to  note  the  difference  between  the  number  of 
modes  required  to  reconstruct  the  flow  and  the  number  of  modes  required  for  effective  low¬ 
dimensional  modeling  for  control  design.  For  real-time  control,  we  are  interested  in  estimating  only 
those  modes  required  for  closed-loop  control.  On  the  other  hand,  an  accurate  reconstruction  of  the 
velocity  field  based  on  a  low-dimensional  model  may  be  obtained  using  between  4-8  modes. 

A  very  important  lesson  learned  concerns  the  validity  of  the  low-dimensional  model.  Is  there  an 
operating  envelope  wherein  the  model  is  valid?  How  do  we  find  this  envelope?  The  cylinder  wake 
flow  can  be  forced  in  an  open  loop  fashion  using  sinusoidal  displacement  of  the  cylinder  with  a  given 
amplitude  and  frequency.  Koopman  (11967)  investigated  the  response  of  the  flow  to  this  type  of 
forcing  in  a  wind  tunnel  experiment.  He  found  a  region  around  the  natural  vortex  shedding  frequency 
where  he  could^chieve  “lock-in”,  which  is  characterized  by  the  wake  responding  to  the  forcing 
by  establishing  a  fixed  phase  relationship  with  respect  to  the  forcing.  The  frequency  band  around  the 
natural  vortex  shedding  frequency  for  which  lock-in  may  be  achieved  is  amplitude  dependent,  as 
shown  in  Fig.  3.3.  Inside  the  V-shaped  area  the  flow  responds  to  the  forcing  with  a  fixed  phase 
relationship,  outside  the  response  to  the  forcing  is  chaotic. 


3.3:  Validity  of  the  Low-dimensional  Model 


Fig. 
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In  general,  the  larger  the  amplitude,  the  larger  the  frequency  band  for  which  lock-in  is  possible. 
However,  a  threshold  amplitude  exists  below  which  the  flow  will  not  respond  to  the  forcing  any  more. 
In  Koopman’s  experiment,  this  amplitude  was  at  10%  peak  displacement  of  the  cylinder.  Shifting  the 
forcing  frequency  away  from  the  natural  shedding  frequency  yields  a  qualitatively  different  behavior, 
ultimately  yielding  a  chaotic  flow  behavior  at  and  beyond  the  lock-in  limit  according  to  Koopman 
(1967).  We  were  able  to  verify  this  behavior  as  shown  in  Fig.  3.3.  The  open  loop  forcing  results  have 
important  implications  for  the  closed  loop  feedback  control  runs.  Since  our  POD  model  is  based  on 
unforced  flow  field  data,  it  can  only  capture  flow  behavior  that  possesses  similar  phenomenology  when 
compared  to  the  unforced  wake.  In  terms  of  the  lock-in  region,  this  flow  behavior  is  encountered  as 
long  as  the  controller  keeps  the  flow  within  the  lock-in  region.  The  chaotic  behavior  at  off-natural 
frequencies  is  clearly  not  modeled  in  the  POD  modes.  Also,  more  importantly,  if  the  displacement  of 
the  cylinder  becomes  smaller  than  about  5%  of  the  cylinder  diameter,  the  flow  will  no  longer  be 
responsive  to  th?  forcing. 


Time  Coefficients  during  Ramped  Open-Loop  Forcing 


Fig.  3.4:  Time  depending  POD  coefficients  obtained  using  transient  forcing 


In  several  applications  of  the  POD  method  for  low-dimensional  model  building  the  mean  flow 
is  removed.  Practically  speaking,  from  an  ensemble  of  snapshots,  the  'average  snapshot'  is  computed 
and  then  this  profile  is  subtracted  from  each  member  of  the  ensemble.  This  is  done  mainly  for  reasons 
of  scale;  i.e.  the  deviations  from  the  mean  contain  information  of  interest  but  may  be  small  compared 
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with  the  original  signal.  However,  Siegel  et  al  (2003)  demonstrated  that  as  the  closed-loop  controller 
reduces  the  unsteady  vortex  induced  drag  in  the  cylinder  wake,  there  is  a  basic  "recovery"  of  energy  in 
the  mean-flow  mode.  It  is  important  to  note  that  the  desire  of  the  closed-loop  controller  is  not  to 
dissipate  the  unsteady  modes  alone  but  to  recover  the  energy  expended  due  to  the  vortex  shedding. 
Therefore  the  low-dimensional  model  should  include  the  mean  flow  mode.  This  vital  conclusion  has 
also  been  noted  by  Siegel  et  al  (2003)  used  the  estimate  of  this  "dynamic"  mean  flow  mode  to  adapt  the 
phase  of  a  PD  (proportional  differential)  controller  applied  to  the  cylinder  wake  at  Re=100.  By 
ensuring  that  the  controller  was  restricted  within  the  operating  envelope  of  the  POD  model,  this 
approach  resulted  in  a  reduction  of  the  vortex  induced  drag  by  90%. 

Given  the  importance  of  capturing  the  energy  exchange  between  the  mean  flow,  aperiodic 
mode,  and  the  unsteady  periodic  modes,  care  needs  to  be  taken  in  selecting  an  appropriate  snapshot 
ensemble  for  thf  POD  procedure.  If  the  snapshots  are  comprised  from  the  steady  state  regime  alone  as 
done  often  in  recent  years,  then  there  is  an  obvious  difficulty  in  capturing  the  desired  transient  energy 
exchange.  As  an  example,  consider  the  cylinder  wake  problem  at  Re=100.  A  POD  procedure  is  applied 
to  500  “snapshots”  from  steady-state  forcing  at  d/D=0.2  and  ramping  up  until  d/D=0.3.  The  mean  flow 
($1  )  was  included  in  the  POD  procedure.  In  addition  to  the  mean  flow  mode  Oj  and  the  von  Karman 
modes  <D 2  and  O3,  a  shift  mode,  O4  was  also  detected.  The  energy  content  of  the  shift  mode  is  about  a 
fourth  of  the  von  Karman  periodic  modes.  About  99.17%  of  the  entire  kinetic  energy  is  contained 
within  the  first  4  modes.  Once  the  essential  dynamics  has  been  captured,  we  can  go  on  and  obtain  the 
desired  set  of  ordinary  differential  equations,  described  in  Equation  3.  This  procedure  has  been 
described  in  detail  by  Luchtenburg  et  al  (2003). 

3.3  Sensor  Configuration,  Estimation  and  Observability 

The  time  histories  of  the  temporal  coefficients  of  the  POD  model  are  determined  by  applying 
the  least  squares  technique  to  the  spatial  modes  and  the  unforced  flow.  The  intent  of  the  proposed 
strategy  is  that  the  velocity  measurements  provided  by  the  sensors  are  processed  by  the  estimator  to 
provide  the  estimates  of  the  first  two  temporal  modes.  The  estimation  scheme,  based  on  the  linear 
stochastic  estimation  procedure  introduced  by  Adrian  (1977),  predicts  the  temporal  amplitudes  of  the 
first  four  POD  modes  from  a  finite  set  of  velocity  measurements  obtained  from  the  CFD  solution  of  the 
uncontrolled  cylinder  wake.  For  each  sensor  configuration,  velocity  measurements,  equally  spaced  at 
an  appropriate  time  interval,  were  used.  The  mode  amplitudes,  01-04,  were  mapped  onto  the  extracted 
sensor  signals,  us,  as  follows:  ' 

m 

o,(')=£c>.(')  <3-4> 


where  m  is  the  number  of  sensors  and  Cs  represents  the  coefficients  of  the  linear  mapping.  The 
effectiveness  of  a  linear  mapping  between  velocity  measurements  and  POD  states  has  been 
experimentally  validated  by  Siegel  et  al.  (2004)  The  coefficients  Cs  (n=l,2;  s=T,m)  in  Equation  (4) 
are  obtained  via  the  linear  stochastic  estimation  method  from  the  set  of  discrete  sensor  signals  and 
temporal  mode  amplitudes. 
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For  the  sensor  configuration,  the  effectiveness  of  the  linear  stochastic  estimation  process 
for  the  estimation  of  the  first  four  temporal  mode  amplitudes,  aj  -  04,  is  calculated.  The  quantitative 
criterion,  associated  with  estimation  accuracy,  is  based  on  the  root  mean  square  (RMS)  of  the  error 
between  the  estimated  and  the  extracted  mode  amplitudes. 

The  spatial  eigenfunctions  obtained  from  the  POD  procedure  provides  information  concerning 
the  location  of  areas  where  modal  activity  is  at  its  highest.  These  energetic  areas  would  be  the  maxima 
and  minima  of  the  spatial  eigenfunctions  (Cameron  et  al,  2004).  Placing  sensors  at  the  energetic 
maxima  and  minima  of  each  mode  is  the  basic  hypothesis  of  this  effort  and  the  purpose  of  the  CFD 
simulation  is  to  design  a  sensor  configuration  which  is  later  validated  using  water-tunnel  experiments. 
Location  of  vorticity  maxima/minima  of  the  spatial  eigenfunctions  are  used  for  sensor  placement.  A 
sensor  was  placed  on  each  the  maxima  and  the  minima  of  modes  1  and  2.  On  the  other  hand,  for 
effective  estimation,  two  pairs  of  sensors  each  are  located  on  maxima/minima  of  modes  3  and  4.  The 
estimated  versus  desired  mode  amplitude  plot,  for  the  above  sensor  configuration  is  presented  in  Fig. 
3.5.  Furthermore,  the  details  concerning  the  validity  of  the  measurement  equation  and  conditions  for 
observability  are  presented  in  Siegel,  Cohen  and  McLaughlin  (2002).  The  main  focus  of  the  work  to 
date  of  the  sensor  placement  studies  has  been  the  usage  of  a  steady-  state  scenario.  However,  the 
development  of  an  observer  for  real-time  estimation  of  the  coefficients  during  a  closed-loop 
simulation/experiment  remains  a  challenge. 


100 


CFD  Data  Re  =100 


Time  [s] 


3.2 


Fig.  3.5:  Extracted  Mode  Amplitudes  projected  on  the  estimated  ones  for  CFD  data  at  Re=100 
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3.4  Controllability,  and  Stability  (COS)  Algorithms 

Let  us  examine  Equation  (3.3).  For  the  origin,  o*  =  0  as  the  desired  equilibrium  or  fixed  point, 
the  aim  is  to  find  an  appropriate  control  law  for  fa  that  will  make  the  equilibrium  stable.  The  current 
effort  only  concentrates  on  stability  issues.  The  function  gk  can  be  expanded  locally  as  a  Taylor  series 
about  the  desired  equilibrium  point: 


OCtj  (3.5) 

where  °(\ak  Represents  higher  order  terms  of  the  expansion  and  can  be  neglected  for  the  stability 

analysis.  For  the  wake  POD  model,  we  observe  gk(0)  ~  0.  Therefore,  the  linearization  of  ak  about  the 
desired  equilibrium  point  is: 


ak  —  J  -  a  +bK  -  fa 


where  the  Jacobian,  J,  is: 

"  dg,(0)  fo(0) 

dax  daM 

-  Sg'u\0)  %>) 

da,  daM 

For  nonlinear  system  given  in  Equation  (3.3),  the  simplest  approach  to  study  controllability  is 
to  consider  its  linearization  as  described  in  Equation  (3.6). 

Definition:  The  pair  (J,  B)  is  state  controllable  if  and  only  if  there  exists  a  control  fa  that  will  transfer 
any  initial  state  Ak(t  =  0)to  the  desired  equilibrium  point  in  finite  time. 

The  following  algebraic  condition  for  controllability  may  be  written  down  for  B  =  [bj,b2,...bM]: 

rank(B:JB:J2B:-:J1’_1B)  =  n  (37) 


21 

STTR  Final  Report,  July  2004 


Cobalt  Solution,  LLC 
USAF  Academy 
STTR  Topic  AF03T007 
Contract:  F49620-03-C-0108 


t  0  l  it  t  i  Q  »  1 


The  conditions  for  asymptotic  stability  can  now  be  stated  for  linearized  models  about  their 
equilibrium  point,  as  follows:  For  the  linearization  given  in  Equation  (3.6)  if  the  Jacobian,  Jc ,  has  n 
eigenvalues,  each  of  which  has  a  strictly  negative  real  part,  then  the  equilibrium  point  is 
asymptotically  stable. 

In  addition  to  conditions  for  stability,  it  is  also  important  to  make  sure  that  the  closed-loop 
linearized  system  is  hyperbolic.  A  fixed  point  of  an  nth  order  system  is  hyperbolic  if  all  the 
eigenvalues  of  the  linearization  (Jacobian)  lie  off  the  imaginary  axis.  The  Hartman-Grobman  theorem 
states  that  the  local  phase  portrait  near  a  hyperbolic  fixed  point  is  “topologically  equivalent”  to  the 
phase  portrait  of  the  linearization;  in  particular  the  stability  type  of  the  fixed  point  is  faithfully  captured 
by  the  linearization  (Cohen,  2004).  A  linearized  system,  see  Equation  (3.6),  that  is  hyperbolic  is 
equivalent  in  terms  of  stability  and  bifurcations,  chaos  and  attractors,  equilibria  and  limit  cycles  to  the 
nonlinear  POEffnodel  (see  Equation  (3.3)).  From  a  practical  point  of  view,  let  us  consider  the  case 
when  the  controller  is  simply  based  on  proportional  control  feeding  back  on  the  estimate  of  mode  1 
alone,  namely: 


fa  ~  ~KP  ■  a\st  <»> 

where  Kp  is  the  proportional  gain  of  the  P  controller  and  aiest  is  the  estimate  of  the  time-dependent 
coefficient  of  Mode  1,  aj.  Inserting  the  control  law  in  Equation  (3.8)  into  Equation  (3.6)  of  yields: 


dk=Jc-  a, 


(3.9) 


where  Jc  is  the  “closed-loop”  Jacobian  and  a  linear  stability  analysis  based  on  Jc  will  provide  an 
insight  into  the  behavior  of  the  closed-loop  system: 


J  c  — 


Q) 

dA, 

dgA  0) 

dA , 


-bvKf 


-K-Kr 


3g,(  0) 
Mu 

Sgy(0) 


t 


Now,  it  is  the  aim  of  the  control  design  to  find  an  appropriate  gain,  Kp,  which  will  render  all 
the  eigenvalues  of  Jc  to  have  a  negative  real  part.  In  addition,  the  eigenvalues  need  to  lie  off  the 
imaginary  axis  by  an  adequate  margin  so  that  the  system  is  hyperbolic.  An  illustration  of  this  approach 
for  the  cylinder  wake  problem  is  detailed  in  Cohen  et  al  (2003). 
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3.5  Important  Control  Lessons  Learned 

Table  6.1  summarizes  the  important  lessons  learned  by  the  multi-disciplinary  research  team  at 
the  Department  of  Aeronautics  of  the  Air  Force  Academy  in  the  context  of  the  cylinder  wake  problem. 
These  lessons  result  from  the  application  of  a  unique  suite  of  simulation  and  experimental  tools 
development  by  the  team.  The  path  towards  meaningful  applications  of  closed-loop  control  has 
several  obstacles  that  need  to  be  effectively  addressed.  These  open  issues,  also  listed  in  Table  6.1,  will 
be  continued  to  be  addressed  as  we  proceed  towards  Phase  II. 


Description 

Lessons  Learned 

Open  Issues 

* 

Modeling 

•  Validity  of  low-dimensional  model  connected 
to  lock-in  phenomena 

•  Energy  exchange  between  aperiodic  "mean 
flow"  modes  and  periodic  vortex  shedding 
modes 

•  Snap-shot  ensemble  to  include  data  reflecting 
transient  forcing 

•  Evaluate  sensitivity  of  low¬ 
dimensional  model  using 
transient  forcing  for 
Reynolds  number  effects, 
and  geometrical  different 
bluff  bodies 

Sensor 

Placement  and 
Number 

•  Sensor  placement  may  be  based  on  the 
topology  of  the  spatial  Eigenfunctions 

•  Number  of  sensors  required  are  determined  by 
keeping  RMS  error  of  estimated  signal  within  a 
desired  threshold 

•  Modification  of  strategy 
using  transient  forcing 

•  Modify  approach  for 
surface  mounted  sensors 
as  opposed  to  sensors  in 
the  flow 

Observability, 

Controllability 

Conditions 

•  Controllability  based  on  the  linearization  of  the 
state  equation 

•  Observability  is  based  on  linearized  state 
equation  and  the  measurement  equation 

•  Validate  approach  using 
transient  forcing  model 

Estimator 

•  Linear  stochastic  estimator  effectively  provides 
measurement  equation  from  steady-state  data 

•  Modify'' observer  to 
effectively  operate  with 
few  sensors  in  a  transient 
environment 

Controller 

•  Feedback  is  effective  when  based  on  primary 
shedding  mode  of  von  Karm&n  vortex-street 

•  Adaptive  strategy,  based  on  the  estimation  of 
the  "mean-flow"  mode  incorporated  to  tune  the 
phase  of  a  Proportional-Differential  (PD) 
controller 

•  Experimental  validation 
of  adaptive  controller 
that  provides  augmented 
performance  using  the 
.  US  Air  Force  Academy 

cylinder  wake 
benchmark 

Table  6.1:  Control  Issues:  Summary  of  Lessons  Learned  and  Open  Issues 
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4.  Application  I:  Circular  Cylinder  Wake  Control 

The  main  portion  of  this  section  was  accomplished  under  core  AFOSR  funding. 
Nonetheless,  the  results  are  pertinent  to  the  toolbox  development  and  approach  taken  for  all  three 
applications.  In  particular,  the  feedback  control  strategy  used  for  the  circular  cylinder  wake  was 
carried  on  to  the  D  shaped  cylinder  presented  in  the  next  section. 

4.1  Circular  Cylinder:  Introduction 

Two-Dimensional  bluff  body  wakes  have  been  investigated  for  quite  some  time.  In  a  two- 
dimensional  cylinder  wake,  self-excited  oscillations  in  the  form  of  periodic  shedding  of  vortices 
are  observed  above  a  critical  Reynolds  number  of  approximately  47.  This  behavior  is  referred  to 
as  the  von  Karman  Vortex  Street.  According  to  Williamson,  the  regime  of  laminar  vortex 
shedding  extends  to  a  Reynolds  number  of  approximately  180,  before  three-dimensional 
instabilities  occur.  This  is  the  Reynolds  number  regime  that  we  target  in  this  investigation. 
However,  the  Karman  vortex  street  as  the  fundamental  feature  of  this  type  of  wake  flow  is 
sustained  to  very  large  Reynolds  numbers  (on  the  order  of  millions).  Therefore  the  lessons 
learned  at  low  Reynolds  numbers  will  still  be  applicable  to  applications  of  practical  interest  at 
much  higher  Reynolds  numbers.  Conversely,  it  would  be  impossible  to  control  the  flow  at  high 
Reynolds  numbers  without  being  able  to  successfully  do  this  at  low  Reynolds  numbers. 

The  non-linear  oscillations  of  the  vortex  street  lead  to  some  undesirable  effects  associated 
with  unsteady  pressures  such  as  fluid-structure  interactions  and  lift/drag  fluctuations.  Also,  the 
vortices  themselves  greatly  increase  the  drag  of  the  bluff  body,  compared  to  the  steady  wake  that 
can  be  observed  at  lower  Reynolds  numbers.  Monkewitz  showed  that  the  von  Karman  Vortex 
Street  is  the  result  of  an  absolute,  global  instability  in  the  near  wake  of  the  cylinder.  Further 
downstream  the  flow  is  convectively  unstable.  This  absolute  instability  is  causing  the  flow  to 
behave  as  a  self  sustained  oscillator,  with  internal  positive  feedback  leading  to  temporal 
amplification  of  the  oscillation  by  the  recirculation  region  downstream  of  the  cylinder. 

Many  attempts  to  improve  the  unsteady  vortex  street  have  been  made.  When  active  open 
loop  forcing  of  the  wake  is  employed,  the  vortex  shedding  in  the  wake  can  be  locked  in  phase  to 
the  forcing  signal.  While  these  findings  suggest  that  the  dominant  structures  in  the  flow  field  can 
be  influenced  by  the  forcing,  it  also  strengthens  the  vortices,  and,  consequently  increases  the 
mean  drag  as  well  as  unsteady  lift  fluctuations.  Different  forcing  methods  are  effective  in 
influencing  the  behavior  of  the  flow.  Acoustic  excitation  of  the  wake,  longitudinal,  lateral  or 
rotational  vibration  of  the  cylinder  model,  and  alternate  blowing  and  suction  at  the  separation 
points  have  been  used.  Using  these  methods,  the  flow  exhibits  regions  of  lock-in  and  non-lock- 
in.  Koopmann  experimentally  investigated  these  regions,  and  found  that  the  lock-in  frequency 
range  depends  on  the  forcing  amplitude.  The  higher  the  forcing'  amplitude,  the  larger  the 
frequency  band  for  which  he  could  achieve  lock-in.  Additionally,  even  at  the  natural  vortex 
shedding  frequency,  he  found  a  minimum  threshold  amplitude  that  was  needed  for  lock-in  to 
occur. 
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All  of  these  open  loop  forcing  methods  have  not  been  shown  to  reduce  the  drag, 
independent  of  frequency  and  amplitude  employed.  The  only  exceptions  are  situations  where  the 
separation  point  location  is  shifted.  It  should  be  noted  that  the  geometry  of  a  circular  cylinder 
lends  itself  to  active  control  methods  that  target  the  separation  point  location  of  the  flow  rather 
than  the  absolute  instability  of  the  wake  itself.  Using  methods  like  periodic  blowing  and  suction 
on  the  cylinder  surface  in  a  fashion  similar  to  that  employed  on  the  suction  side  of  airfoils,  the 
separation  point  can  be  moved  aft  which  in  turn  will  lead  to  a  narrower  wake.  A  narrower  wake 
will  exhibit  improved  stability  characteristics,  in  addition  to  lower  drag  due  to  a  lower  velocity 
deficit  in  and  by  itself.  This  effect  can  be  observed  in  the  natural  cylinder  wake  during  the  “drag 
crisis”,  when  laminar-turbulent  boundary  layer  transition  occurs  upstream  of  the  separation 
point,  and  the  resulting  turbulent  boundary  layer  keeps  the  flow  attached  longer  shifting  the 
separation  point  downstream.  Thus  feedback  control  investigations  using  periodic  blowing  and 
suction  like^hat  employed  by  Min  and  Choi  actually  employ  two  flow  control  techniques 
simultaneously,  namely  separation  control  and  wake  stabilization  due  to  feedback.  It  is  difficult 
if  not  impossible  to  judge  what  portion  of  the  improvement  is  due  to  either  of  these  techniques  in 
their  simulations. 

The  only  way  of  suppressing  the  self-excited  flow  oscillations  without  altering  the  mean 
flow  is  by  the  incorporation  of  active  closed-loop  flow  control.  Traditionally,  several 
fundamentally  different  approaches  to  achieve  feedback  flow  control  have  been  employed.  The 
mathematically  driven  approach  to  develop  a  control  scheme  is  hampered  by  the  complexity  of 
the  governing  Navier  Stokes  equations.  In  order  to  tackle  this  complexity,  one  needs  to  make 
simplifying  assumptions.  At  this  point,  the  assumptions  made  in  simplifying  the  equations  have 
often  rendered  the  results  inapplicable  to  real  life  experiments  (Li  and  Aubry).  If  no 
simplifications  or  assumptions  are  made,  however,  the  resulting  control  algorithm  (if  it  can  be 
derived  at  all  with  today’s  computing  capabilities)  is  often  too  complex  to  be  implemented  in 
real  time  (Bewley  and  Trenchea). 

On  the  other  hand,  approaching  the  controls  problem  using  an  experimental  /  empirical 
approach  without  any  modeling  of  the  physics  of  the  flow  yields  unsatisfactory  results  also. 
Experimental  studies  conducted  by  Roussopulos  show  that  a  linear  proportional  feedback  control 
based  on  a  single  sensor  feedback  is  able  to  delay  the  onset  of  the  wake  instability  only  slightly, 
rendering  the  wake  stable  at  Re  about  20%  higher  than  the  unforced  case.  Roussopulos  reports 
that  above  Re  =  60,  a  single-sensor  feedback  may  suppress  the  original  mode  but  destabilizes  one 
of  the  other  modes.  Therefore,  better  control  strategies  are  needed  to  stabilize  the  wake  at 
Reynolds  numbers  of  technical  interest. 

The  solution  to  this  problem  lies  in  the  development  of  a  low  order  model  of  the  flow. 
The  model  can  be  used  both  for  controller  development,  as  well  as  flow  field  state  estimation. 
Ideally,  it  reduces  the  complexity  of  the  governing  Navier  Stokes  equations  to  a  level  that  the 
model  can  be  implemented  in  real  time,  while  still  capturing  the  important  physics  of  the  flow. 
Gillies  pioneered  the  application  of  this  technique  to  cylinder  wakes  by  developing  a  low 
dimensional  model  of  the  cylinder  wake  at  a  Reynolds  number  of  100.  Cohen  et  al.  have  shown 
that  using  this  model,  the  cylinder  wake  model  flow  can  be  successfully  controlled  using  a 
relatively  simple  linear  control  approach  based  on  the  most  dominant  mode  only. 
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The  goal  of  this  effort  is  to  apply  the  approach  developed  by  Cohen  et  al.  to  a  full  Navier 
Stokes  simulation  of  the  flow  field.  A  limited  number  of  sensors  placed  in  the  wake  are  used  to 
estimate  the  state  of  the  flow  which  is  characterized  using  a  low  dimensional  model.  The 
controller  then  acts  on  the  flow  state  estimates  in  order  to  determine  the  actuator  displacement 
(Figure  4.1  shows  the  overall  setup  of  this  experiment). 


Figure  4.1  Feedback  Control  Setup 


4.2  Circular  Cylinder:  Numerical  Methods 
CFD  Model 

For  the  numerical  simulations,  Cobalt  Solutions’  Cobalt  solver  V.2.02  was  used.  This 
code  can  operate  in  many  different  modes  using  various  turbulence  models.  However,  for  the 
present  investigation  it  was  used  as  a  direct  Navier  Stokes  solver  with  second  order  accuracy  in 
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time  and  space.  Cobalt  operates  on  unstructured  grids.  For  all  investigations  an  unstructured  two- 
dimensional  grid  with  63700  nodes  /  31752  elements  was  used,  see  Figure  4.2.  The  grid  extended 
from  -16.9  cylinder  diameters  to  21.1  cylinder  diameters  in  the  x  (streamwise)  direction,  and 
±19.4  cylinder  diameters  in  y  (flow  normal)  direction. 

Other  pertinent  simulation  parameters: 

•  Reynolds  Number  (Re)  =100 

•  Mean  flow,  U=  34  m/s 

•  Damping  Coefficients:  Advection  =  0.01;  Diffusion  =  0.00 

•  32  Iterations  for  matrix  solution  scheme 

•  3  Newtonian  sub-iterations 

•  Non-dimensional  time  step  At*  =  At*U/D=0.05 

•  Time  step,  At  =  0.00147  s 

•  Vortex  shedding  frequency  5.55  Hz 


The  simulation  was  triggered  by  skewing  the  incoming  mean  flow  by  a  =  0.5  degrees  in 
order  to  introduce  an  initial  perturbation.  An  important  issue  concerning  the  validity  of  the  CFD 
model  needs  to  be  addressed  before  using  the  data.  A  grid  and  time  resolution  study  showed 
good  convergence  for  the  simulation  parameters  outlined  above.  For  further  validation  of  the 
unforced  cylinder  wake  CFD  model  at  Re  =  100,  the  resulting  value  of  the  mean  drag  coefficient, 
cd,  was  compared  to  experimental  and  computational  investigations  reported  in  the  literature.  At 
Re  =  100,  experimental  data,  reported  by  Oertel  and  Panton  point  to  c,  values  ranging  from  1.26 
to  1.4.  Furthermore,  Min  and  Choi  report  on  several  numerical  studies  that  obtained  drag 
coefficients  of  1.35  and  1.337.  The  COBALT  CFD  model,  used  in  this  effort  results  in  a  cd  =1.35 
at  Re  =  100,  which  compares  well  with  the  reported  literature.  Another  important  benchmark 
parameter  is  the  non-dimensional  vortex  shedding  frequency,  the  Strouhal  number  (St)  for  the 
unforced  cylinder  wake.  Experimental  results  at  Re  =  100,  presented  by  Williamson,  show 
Strouhal  numbers  ranging  from  0.163  to  0.166.  The  Strouhal  number  obtained  from  the 
COBALT  CFD  model  used  in  this  effort  is  St  =  0.163  at  Re  =  100,  which  also  compares  well. 

The  simulations  were  performed  on  a  Beowulf  Linux  cluster  consisting  of  64  Pentium  III 
processors  operating  at  1  Ghz.  When  running  on  8  processors,  typically  a  time  step  took  on  the 
order  of  6  s  to  compute.  Employing  larger  number  of  processors  yielded  disproportionately  small 
improvements  in  execution  time,  due  to  network  and  disk  access  overhead  for  saving  the  results 
at  the  end  of  each  time  step. 

4.3  Circular  Cylinder:  POD  Modeling  and  Estimation 

Feasible  real  time  estimation  and  control  of  the  cylinder  wake  may  be  effectively  realized 
by  reducing  the  model  complexity  of  the  cylinder  wake  as  described  by  the  Navier-Stokes 
equations,  using  POD  techniques.  POD,  a  non-linear  model  reduction'approach  is  also  referred  to 
in  the  literature  as  the  Karhunen-Loeve  expansion.  The  desired  POD  model  contains  an  adequate 
number  of  modes  to  enable  modeling  of  the  temporal  and  spatial  characteristics  of  the  large-scale 
coherent  structures  inherent  in  the  flow. 
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In  this  effort,  the  method  of  “snapshots”  introduced  by  Sirovich  is  employed  to  generate 
the  basis  functions  of  the  POD  spatial  modes  from  the  numerical  solution  of  the  Navier-Stokes 
equations  obtained  using  COBALT.  In  all  200  snapshots  were  used,  equally  spaced  0.00735  s 
apart.  The  time  between  snapshots  is  five  times  the  simulation  time  step.  The  snap-shots  were 
taken  after  ensuring  that  the  cylinder  wake  reached  steady  state.  Only  the  U  velocity  component 
(in  the  direction  of  the  mean  flow)  was  used  for  POD  decomposition  in  this  effort.  This  decision 
was  made  in  order  to  be  able  to  estimate  the  mode  amplitudes  based  on  sensor  information, 
which  in  future  experiments  will  yield  the  U  and  V  component  of  velocity.  Since  the  change  in 
mean  flow  distribution  is  an  important  quantity,  we  chose  the  U  velocity  component  over  the  V 
velocity  component.  We  found  that  more  than  99.98%  of  the  kinetic  energy  of  the  flow  lies  in 
the  first  eight  modes,  with  more  than  93.5%  in  the  first  four  modes.  An  important  aspect  of 
reduced  order  modeling  concerns  truncation:  how  many  modes  are  important  and  what  are  the 
criteria  for  effective  truncation? 

The  answers  to  the  above  questions  have  been  addressed  by  Cohen  et  al.  This  effort 
showed  that  control  of  the  POD  model,  of  the  von  Karman  vortex  street  in  the  wake  of  a  circular 
cylinder  at  Re  =  100,  is  enabled  using  just  the  first  mode.  Furthermore,  feedback  based  on  the 
first  mode  alone  suppressed  all  the  other  modes  in  the  four  mode  POD  model,  indicating  that 
higher  order  modes  derive  from  the  fundamental  modes.  In  view  of  this  result,  truncation  of  the 
POD  model  took  place  after  the  first  four  modes.  At  this  point,  it  is  imperative  to  note  the 
difference  between  the  number  of  modes  required  to  reconstruct  the  flow  and  the  number  of 
modes  required  to  control  the  flow.  In  this  effort,  we  are  interested  in  estimating  only  those 
modes  required  for  closed-loop  control.  On  the  other  hand,  an  accurate  reconstruction  of  the 
velocity  field  based  on  a  low-dimensional  model  may  be  obtained  using  between  4-8  modes.  The 
POD  algorithm  was  applied  to  the  fluctuating  velocity  component  in  the  direction  of  the  flow  (u) 
as  described  in  Equation  (1).  The  decomposition  of  this  component  of  the  velocity  field  is  as 
follows: 


u(x,  y,  t)  -  U(x,  y)  +  u(x,  y,  t) 


where  U[m/s]  denotes  the  mean  flow  velocity  and  u[m/s]  is  the  fluctuating  component  that  may 
be  expanded  as: 


u(x,  y,  t)  =  f]  ak  (t)<J)[k)  (x,  y) 

k=1  (4.2) 

where  al(t)  denotes  the  time-dependent  coefficients  having  units  of  m/s  and  <j)(x,  y)  represent  the 
non-dimensional  spatial  Eigenfunctions  (see  Figure  4.3)  determined  from  the  POD  procedure. 
Shown  are  the  first  four  modes  of  the  POD  decomposition,  plus  a  5th  mode  that  was  derived  by 
subtracting  the  mean  ffeestream  velocity  from  the  mean  flow  distribution  of  the  unforced  flow 
field.  This  mode  was  found  to  be  necessary  to  obtain  an  estimate  of  the  effect  of  feedback  flow 
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Figure  4.3  Eigenfunctions  of  the  5  Mode  POD  model  based  on  CFD  data  at  Re  =  100,  using  the 
U  velocity  component  as  input  for  the  POD  decomposition.  Flow  from  left  to  right,  the  cylinder 
is  located  at  (0,0)  and  has  a  diameter  of  1 


Once  the  spatial  POD  eigenfunctions  have  been  derived,  the  corresponding  time- 
dependent  coefficients  al(t),  or  mode  amplitudes,  need  to  be  calculated.  For  this,  two  different 
schemes  are  reported  in  literature.  Most  often  a  Galerkin  projection  is  used,  which  involves 
projecting  the  spatial  eigenfunctions  onto  the  Navier  Stokes  equati'bns.  This  process  involves 
spatial  derivatives  of  the  snapshots,  which  are,  particularly  in  the  case  of  experimental  data, 
inherently  sensitive  to  noise.  Gillies8  used  a  simple  least  squares  fit,  which  we  found  to  be  much 
more  robust.  While  we  employ  the  least  square  fit  in  the  CFD  simulations,  the  experiment  will 
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make  use  of  linear  stochastic  estimation  (LSE)  in  order  to  estimate  the  mode  amplitudes  in  real 
time.  LSE  is  deterministic  in  terms  of  computing  time,  while  least  square  fitting  is  not.  Thus  LSE 
is  a  much  better  choice  for  real  time  implementation. 

For  the  feedback  controlled  runs,  the  CFD  solver  writes  sensor  information  at  requested 
(x,  y)  locations  in  the  flow  to  a  file  after  calculating  a  time  step,  and  then  waits  for  an  external 
control  algorithm  to  update  the  file  with  the  new  control  command  for  the  next  time  step.  The 
sensor  grid  employed  for  all  simulations  is  shown  in  Figure  4.2  and  employs  a  total  of  35  sensors 
in  the  near  wake  of  the  cylinder.  The  main  advantage  of  this  sensor  grid  over  others  investigated 
is  in  its  ability  to  provide  a  global  estimate  of  the  mode  amplitudes  that  shows  little  error 
compared  to  using  all  grid  points.  This  holds  true  both  for  the  unforced  case  as  well  as  the 
feedback  controlled  cases.  Typical  errors  are  negligible  in  phase  and  less  than  5  %  in  amplitude. 

Whirf  no  extensive  efforts  to  optimize  the  sensor  locations  was  undertaken,  we  compared 
a  localized  uniformly  spaced  sensor  field  with  35  sensors  between  x/D  =  0.75  and  x/D  =  1.75  to 
the  more  widely  spaced  confiyguration  shown  in  Figure  4.2.  While  both  grids  performed  equally 
well  in  estimating  the  unforced  flow  field,  for  the  feedback  control  runs  the  local  grid  develops 
large  phase  and  amplitude  errors  as  soon  as  the  flow  responds  to  the  forcing.  With  this  finding 
we  decided  to  use  the  distributed  sensor  field  shown  in  Figure  4.2. 
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Figure  4.2  Computational  grid  (dots)  and  sensor  locations  (circles). 
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4.4  Circular  Cylinder:  Controller 

Figure  1  shows  the  controller  block  diagram.  The  main  components  are  a  Mode 
estimator,  a  Proportional  and  Differential  (PD)  controller  and  a  low  pass  filter  with  a  comer 
frequency  of  four  times  the  natural  shedding  frequency.  The  Cobalt  CFD  solver  has  the  ability  to 
perform  rigid  body  motion  of  a  given  grid.  This  feature  was  used  to  perform  both  time 
periodically  forced  and  feedback  controlled  simulations  with  one  degree  of  freedom.  For  all 
investigations,  only  displacement  of  the  cylinder  in  flow  normal  (y)  direction  was  employed  for 
forcing  the  flow.  The  control  algorithm  acts  on  the  estimate  of  the  Mode  1  amplitude  only.  This 
design  decision  was  made  based  on  our  earlier  investigations  controlling  a  low  dimensional 
model  of  the  flow.  For  the  low  dimensional  model,  proportional  gain  applied  to  Mode  1  only  was 
sufficient  to  suppress  vortex  shedding.  Since  our  CFD  simulations  require  a  filter  to  avoid 
feeding  backfsmall  amounts  of  noise,  we  employed  a  PD  feedback  control  strategy: 


y^-K^  +  K, 


dax 


dt  (4.3) 

Instead  of  directly  specifying  the  Kpand  Kd  gains,  these  can  be  expressed  in  terms  of  an 
overall  gain  K  and  a  phase  advance  cp: 


K  =K-cos(<p) 

K  = 

2  nf  (4.4) 

With  f  being  the  natural  vortex  shedding  frequency. 

Physically,  the  control  algorithm  was  implemented  in  Matlab  on  a  separate  PC  running 
Windows.  It  interfaced  to  the  Beowulf  cluster  running  Cobalt  using  Windows  file  sharing 
through  Samba,  to  read  the  sensor  information  and  update  the  cylinder  displacement. 

4.5  Circular  Cylinder:  Results 

Before  closed  loop  feedback  flow  control  is  employed,  it  is  important  to/ investigate  the 
dynamics  of  the  unforced  flow  field  in  detail.  Equally  important,  the  effect  of  open  loop  forcing 
needs  to  be  understood,  since  the  receptivity  of  the  flow  to  forcing  will  manifest  itself  in  these 
investigations.  The  following  section  will  outline  the  results  of  these  investigations,  and  also 
show  the  limitations  of  the  type  of  forcing  employed  as  well  as  the  limits,  of  the  flow 
improvements  that  may  be  obtained  using  feedback  control. 

The  following  two  sections  will  highlight  a  few  select  cases  of  the  unforced  flow,  open 
loop  forced  flow  and  feedback  controlled  flow.  Two  different  kinds  of  feedback  control  were 
employed,  one  using  a  fixed  set  of  proportional  and  differential  gains,  and  one  set  where  the 
gains  were  varied  depending  on  the  change  in  the  mean  flow.  ‘The  latter  is  usually  referred  to  as  a 
gain  scheduling  scheme. 


31 

STTR  Final  Report,  July  2004 


Cobalt  Solution,  LLC 
USAF  Academy 
STTR  Topic  AF03T007 
Contract:  F49620-03-C-0108 

Unforced  Wake  Properties 

In  a  CFD  simulation,  the  flow  field  is  started  abruptly  at  time  zero.  Figure  4.4  shows  the 
evolution  of  the  flow  field  after  this  startup.  The  flow  evolves  from  a  Stokes-type  streamline 
pattern  at  t  =  0  s  (Figure  4.4  (a))  through  a  steady  wake  with  two  closed  recirculation  bubbles  at  t 
=  0.4  s  (Figure  4.4  (b))  into  the  steady  state  showing  the  unsteady  von  Karmdn  Vortex  Street  at  t 
=  2.94  s  (Figure  4.4  (d)).  During  this  startup,  the  flow  reaches  a  state  of  minimum  drag  at  t  =0.7  s 
(Figure  4.4  (c)).  This  minimum  drag  coincides  with  a  maximum  amplitude  of  mode  5  (the  free 
stream  mode)  as  well  as  a  maximum  mean  recirculation  zone  length  with  the  downstream  end  of 
the  recirculation  zone  located  at  x/D  =  5.4  (not  depicted).  It  is  worth  noting  that  the  minimum 
drag  does  not  coincide  in  time  with  the  steady  wake  as  one  might  expect,  but  rather  with  a  vortex 
shedding  pattern  with  a  very  large  wavelength,  as  shown  in  Figure  4.4  (c).  The  total  drag  in  this 
situation  is  «t>out  16  %  less  than  in  the  steady  state  vortex  shedding  situation.  Thus  one  may 
argue  that  a  feedback  control  scheme  aiming  to  suppress  the  vortex  shedding  may  be  able  to 
recover  up  to  this  portion  of  the  total  drag,  at  best.  We  refer  to  this  16%  portion  of  the  overall 
drag  force  as  the  vortex  induced  drag,  since  it  is  caused  by  the  vortex  shedding  in  the  unsteady 
wake  flow.  It  is  a  portion  of  the  pressure  drag.  About  2  seconds  after  the  startup  of  the 
simulation,  the  wake  approaches  a  time  periodic  vortex  shedding  state.  The  mean  recirculation 
zone  ends  at  x/D  =  1 .9  in  this  final  flow  state. 
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Figure  4.4  (a)-(d)  Instantaneous  streamlines  at  t  =  0.0,  0.4,  0.7  and  3.0  seconds  after  startup  of 
the  simulation,  (e)  lift  and  drag  during  startup  of  the  simulation,  (f)  Mode  amplitudes  during 
startup  of  the  simulation,  for  spatial  mode  distributions  see  Figure  3. 


Open  Loop  Forced  Wake  Properties 

The  cylinder  wake  flow  can  be  forced  in  an  open  loop  fashion  using  sinusoidal 
displacement  of  the  cylinder  with  a  given  amplitude  and  frequency.  Koopman  investigated  the 
response  of  the  flow  to  this  type  of  forcing  in  a  wind  tunnel  experiment.  He  found  a  region 
around  the  natural  vortex  shedding  frequency  where  he  could  achieve  “lock-in”,  which  is 
characterized  by  the  wake  responding  to  the  forcing  by  establishing  a  fixed  phase  relationship 
with  respect  to  the  forcing.  The  frequency  band  around  the  natural  vortex  shedding  frequency  fo 
for  which  lock-in  may  be  achieved  is  amplitude  dependent,  as  shown  in  Figure  4.5(a)  Inside  the 
V-shaped  area  the  flow  responds  to  the  forcing  with  a  fixed  phase  relationship,  with  the  vortex 
shedding  in  phase  with  the  forcing.  Outside,  the  response  to  forcing  is  chaotic.  In  general,  the 
larger  the  amplitude,  the  larger  the  frequency  band  for  which  lock-in  is  possible.  However,  a 
threshold  amplitude  exists  below  which  the  flow  will  not  respond  to  the  forcing  any  more.  In 
Koopman’s  experiment,  this  amplitude  was  at  10%  peak  displacement  of  the  cylinder. 

We  resampled  the  lock-in  region  in  the  CFD  simulation  at  select  amplitude  and  frequency 
pairs  shown  in  Figure  4.5(a).  The  simulations  activated  the  forcing  always  at  the /same  time,  3.3 
seconds  after  the  start  of  the  simulation,  which  resulted  in  the  forcing  being  1 80  degrees  out  of 
phase  with  the  vortex  shedding.  A  typical  run  resulting  in  lock-in  is  depicted  in  Figure  4.5(c)  and 
4.5(d).  Figure  4.5(c)  shows  that  the  flow  field  goes  through  a  transient  phase  before  lock-in  is 
achieved  after  about  9  shedding  cycles.  After  the  transient,  a  fixed  phase  relationship  between 
forcing  and  the  unsteady  lift  force  induced  by  the  vortex  shedding  is  established,  as  shown  in 
Figure  4.5(d).  We  refer  to  the  time  during  which  the  flow  adjusts  to  the  forcing  as  the  settling 
time.  For  comparison.  Figure  4.5(e)  and  4.5(f)  show  the  chaotic  response  of  the  flow  to  a  forcing 
input  characterized  by  a  frequency  /  amplitude  combination  just  outside  the  lock-in  range.  The 
flow  does  not  establish  a  fixed  phase  relationship  to  the  forcing  in  that  case.  A  scan  through 
different  forcing  amplitudes  was  performed  at  the  natural  shedding  frequency  (f/fo  =  1)  with 
amplitudes  ranging  from  1  to  30  %  of  the  cylinder  diameter.  The  settling  times  observed  in  these 
cases  are  shown  in  Figure  4.5(b).  While  the  settling  times  are  roughly  constant  for  forcing 
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amplitudes  between  5  and  10%,  for  smaller  amplitudes  a  drastic  increase  in  settling  times  can  be 
observed.  This  manifests  the  behavior  observed  by  Koopman  around  10%  forcing  amplitude, 
albeit  shifted  towards  somewhat  smaller  amplitudes.  There  are  two  possible  explanations  for  this. 
Koopman  used  spanwise  coherence  as  an  indicator  for  lock-in,  which  may  occur  at  larger 
amplitudes  than  the  local  lock-in  observed  in  our  two-dimensional  simulations.  Additionally,  his 
experiment  was  conducted  in  a  wind  tunnel  environment  which  features  more  mean  flow 
turbulence  than  the  CFD  simulations.  This  would  also  tend  to  increase  the  amount  of  forcing 
needed  to  overcome  the  turbulence  and  achieve  lock-in. 

The  open  loop  forcing  results  have  important  implications  for  the  closed  loop  feedback 
control  runs.  Since  our  POD  model  is  based  on  unforced  flow  field  data,  it  can  only  capture  flow 
behavior  that  is  phenomenologically  similar  to  the  unforced  wake.  In  terms  of  the  lock-in  region, 
this  flow  behavior  is  encountered  as  long  as  the  controller  keeps  the  flow  within  the  lock-in 
region.  The  chaotic  behavior  at  off-natural  frequencies  is  clearly  not  modeled  in  the  POD  modes. 
Also,  more  importantly,  if  the  displacement  of  the  cylinder  becomes  smaller  than  about  5%  of 
the  cylinder  diameter,  the  flow  will  no  longer  be  responsive  to  the  forcing. 
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FIGURE  4.5  (a)  Lock-in  region  adapted  from  Koopman  (1966)  (b)  Settling  time  for  different 
forcing  amplitudes  using  sinusoidal  forcing  at  the  natural  shedding  frequency,  (c-d)  Sinusoidal 
forcing  at  the  natural  frequency  using  an  amplitude  of  5%  of  the  cylinder  diameter  (c)  Lift  and 
Drag  (d)  Phase  plot  no  dimensional  lift  force  vs.  cylinder  motion,  (e-f)  Sinusoidal  forcing  at  75% 
of  the  natural  frequency  using  an  amplitude  of  30%  of  the  cylinder  diameter  (e)  Lift  and  Drag  (f) 
Phase  plot  lift  vs.  cylinder  motion. 


Fixed  Phase  Feedback 

A  series  of  simulations  with  different  phase  advance  cp  were  conducted.  The  overall  gain 
K  was  kept  constant  at  a  value  that  resulted  in  initial  displacements  of  less  than  15%  of  the 
cylinder  diameter,  in  order  to  avoid  strongly  nonlinear  effects  that  have  been  reported  in 
literature  at  higher  amplitudes.  Qualitatively,  the  runs  resulted  in  either  an  increase  or  decrease  of 
the  mode  amplitudes.  A  case  with  a  decrease  in  mode  amplitudes  is  shown  in  Figure  4.6.  During 
the  control  run,  the  global  mode  amplitude  of  Mode  1  decreases  from  a  peak  value  of  around  100 
to  a  peak  value  of  less  than  40.  This  decrease  in  mode  amplitude  does  not  just  apply  to  Mode  1, 
but  also  to  Mode  2  and  the  two  higher  order  modes.  These  findings  are  consistent  with  our 
experience  in  controlling  a  low  order  model  of  the  flow12,  where  a  similar  couplihg  between  the 
modes  could  be  observed.  Considering  the  behavior  of  the  unsteady  lift  and  the  drag,  a  reduction 
in  drag  of  about  14%  of  the  unforced  drag  was  observed,  while  the  unsteady  lift  was  reduced  by 
about  50%.  Comparing  the  drag  reduction  of  14%  to  the  minimum  drag  encountered  during 
startup  of  the  CFD  simulation,  which  is  16%  less  than  the  drag  of  the  unforced  flow  field,  we 
find  that  the  simple  feedback  controller  with  a  fixed  gain  and  phase  recovers  more  than  87%  of 
the  steady  wake  flow  drag  reduction.  However,  all  of  the  fixed  phase  simulations  that  led  to  a 
drag  reduction  were  not  able  to  stabilize  the  flow  at  a  low  drag  value.  In  Figure  4.6(a)  at  a  time 
of  7.5  s,  a  reduction  followed  by  an  increase  in  mode  amplitude  can  be  observed.  This  is  the  first 
indication  of  the  onset  of  an  instability  that  will  ultimately  render  the  flow  field  chaotic,  if  the 
simulation  is  continued.  Inspecting  the  cylinder  displacement  during  this  simulation  (Figure 
4.6(b)),  the  onset  of  the  instability  coincides  with  a  cylinder  displacement  that  has  just  dropped 
below  5%  of  the  cylinder  diameter.  Thus  the  onset  of  the  instability  coincides  with  the  loss  of 
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flow  response  found  in  the  periodically  open  loop  forced  runs  discussed  in  the  previous  section. 
The  instantaneous  vorticity  contours  shown  in  Figure  4.6(c)  show  a  wake  flow  where  the 
vortices  form  further  downstream  of  the  cylinder,  compared  to  the  streamlines  of  the  unforced 
flow  shown  in  Figure  4.4(d).  In  the  unforced  flow,  the  vortices  roll  up  between  1  and  2  diameters 
downstream  of  the  cylinder,  while  in  the  feedback  controlled  situation  the  rollup  occurs  between 
3  and  4  diameters  downstream.  As  a  result  the  reverse  flow  region  is  lengthened,  from  x/D  =1.9 
in  the  unforced  case  to  x/D  =  4.3.  Simultaneously  with  the  lengthening  of  the  recirculation  zone 
we  observe  a  reduction  in  the  vortex  shedding  frequency.  The  correlation  between  the  shedding 
frequency  and  the  reduction  in  drag  is  in  qualitative  agreement  with  the  wake  model  proposed  by 
Ahlborn  et  al. 

A  summary  of  the  effect  of  different  amounts  of  fixed  phase  advance  between  mode  1 
and  the  cylinder  motion  on  both  the  mode  amplitudes  and  the  drag  force  is  shown  in  Figure  4.7. 
While  it  is  apparent  that  the  largest  drag  and  mode  amplitude  reductions  are  achieved  for  a  phase 
advance  of  about  35  degrees,  this  phase  advance  is  also  unstable  over  longer  time  periods.  The 
good  correlation  between  the  mode  amplitude  reduction  and  the  drag  reduction  suggests  a  strong 
link  between  these  quantities.  Also,  it  can  be  seen  that  all  mode  amplitudes  experience  a 
proportional  reduction  which  shows  the  coupling  between  the  modes  to  be  strong.  It  is  also 
interesting  to  note  that  a  zero  degree  phase  advance  has  no  impact  on  the  drag. 


36 

STTR  Final  Report,  July  2004 


i  o  !  *  t  f  o  r>  t 


Cobalt  Solution,  LLC 
USAF  Academy 
STTR  Topic  AF03T007 
Contract:  F49620-03-C-0108 


Figure  4.6  Linear  feedback  of  Mode  1  with  30°  Phase  advance.  The  controller  is  activated  at 
3.03  s  and  deactivated  at  8  s.  (a)  Mode  Amplitudes  (b)  Cylinder  Displacement  and  Frequency  (c) 
Instantaneous  Vorticity  Contours  at  t  =  6.5  s 

With  these  findings  on  the  impact  of  fixed  phase  feedback  on  the  wake,  an  important 
question  to  be  asked  is  if  the  wake  can  be  stabilized  at  a  low  mode  amplitude,  and  if  so,  how?  An 
obvious  parameter  to  adjust  in  order  to  achieve  this  is  the  gain  of  the  controller.  Increasing  the 
gain  leads  to  larger  cylinder  displacements,  which  would  assure  the  ability  to  control  the  wake  at 
small  mode  amplitudes.  However,  we  found  that  larger  gains  do  not  stabilize  the  wake  but  rather 
lead  to  a  low  frequency  oscillation  and  instability  even  when  the  cylinder  displacements  are  kept 
above  5%  of  the  cylinder  diameter.  The  other  parameter  to  be  considered  for  adjustment  is  the 
phase  of  the  feedback  which  we  report  on  in  the  following  section. 


Figure  4.7  Mode  amplitudes  and  drag  force  for  various  phase  advance  angles. 
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Variable  Phase  Feedback 

During  an  investigation  into  different  sensor  configurations,  we  used  a  sensor  field  with 
35  sensors  localized  between  x/D  =  0.75  and  x/D  1.75.  As  was  later  discovered,  this  sensor  field 
developed  a  large  estimation  error  with  respect  to  the  phase  error  of  the  Mode  1  estimate,  relative 
to  an  estimate  based  on  the  entire  flow  field,  as  shown  in  Figure  4.8  (f).  Despite  this  phase  error 
the  wake  stabilized  at  an  overall  drag  reduction  of  about  15%  with  an  unsteady  lift  amplitude 
reduction  of  90%  (Figure  4.8(c)),  compared  to  the  unforced  flow  field.  Inspecting  the  phase 
error,  one  can  see  that  due  to  the  effects  of  the  local  sensor  field  the  phase  advance  is  reduced  to 
almost  zero  in  the  steady  state  case.  This  phase  advance  angle  stabilizes  the  flow  field  at  a  low 
level  of  vortex  shedding,  with  the  recirculation  length  extended  to  x/D  =  3.95,  or  more  than  twice 
the  unforced  length.  The  phase  plot  in  Figure  4.8(d)  demonstrates  that  the  wake  develops  a  new 
limit  cycle  af  reduced  mode  1  amplitude  in  the  steady  state  stabilized  case.  The  vorticity  plot  in 
Figure  4.8(e)  is  very  similar  to  the  fixed  phase  feedback  vorticity  distribution  shown  in  Figure 
4.6(c).  The  recirculation  zone  is  clearly  lengthened  in  the  stabilized  flow,  which  leads  to  a 
decrease  in  the  vortex  shedding  frequency  by  23%  (Figure  4.8(b)). 

While  in  the  run  shown  in  Figure  4.8  the  phase  advance  was  a  result  of  the  dense 
localized  sensor  placement,  the  same  effect  can  be  achieved  using  a  global  sensor  field  like  the 
one  shown  in  Figure  4.2  in  combination  with  a  variable  phase  advance  based  on  the  non  - 
fluctuating  mode  5.  Thus  we  find  that  a  variable  gain  strategy  that  adjusts  the  feedback  gains 
according  to  the  change  in  the  mean  flow  achieves  better  results  than  a  fixed  gain  control 
approach.  The  success  of  the  variable  gain  strategy  is  a  logical  result  given  that  the  initial 
feedback  gain  and  phase  are  based  on  the  non-control  led  flow  field  and  recirculation  length. 
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Figure  4.8  Linear  feedback  of  Mode  1  with  variable  Phase  advance,  (a)  Mode  Amplitudes  (b) 
Cylinder  Displacement  and  Frequency  (c)  Lift  and  Drag  (d)  Phase  between  Cylinder  Position 
and  Mode  1  (e)Instantaneous  Vorticity  Contours  at  t=7.5  s  (f)  Phase  advance  during  the  run 
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The  drag  and  unsteady  lift  force  reduction  manifests  itself  in  a  change  in  the  mean  flow, 
as  well  as  the  RMS  distribution.  Figure  4.9  compares  the  unforced  mean  flow  and  RMS 
distributions  to  those  encountered  in  the  stabilized  state,  between  6  and  8  seconds,  in  the 
feedback  controlled  run.  The  recirculation  zone  length  has  almost  doubled  in  length,  and  the 
peak  in  the  RMS  distribution  is  shifted  from  x/D=2.5  to  x/D  =  5.  Also,  it  can  be  seen  that  the 
controlled  wake  up  to  3  diameters  downstream  of  the  cylinder  is  entirely  steady. 

When  applying  feedback  control,  significant  changes  in  the  mean  flow  field  occur,  as 
shown  in  the  previous  section.  It  is  therefore  of  interest  to  investigate  how  the  stability 
characteristics  of  the  mean  flow  are  modified  as  a  result  of  the  mean  flow  changes.  Linear 
stability  analysis  based  on  numerical  solution  of  the  Orr-Sommerfeld  equations  using  spectral 
methods  (Trefelthen)  was  used  to  analyze  these  changes.  Figure  4.10  compares  the  maximum 
growth  rate  ^f  the  unforced  flow  field  at  a  Reynolds  number  of  100  to  the  steady  state  feedback 
controlled  flow  field  (Run  122).  Despite  the  fact  that  the  near  wake  fluctuations  are  suppressed 
by  the  feedback  as  shown  in  the  previous  section,  the  flow  field  has  become  more  unstable 
beyond  two  diameters  downstream  of  the  cylinder.  Comparing  the  unforced  flow  to  a  stable  flow 
field  at  a  subcritical  Reynolds  number  of  40,  one  can  see  that  the  Karman  vortex  street  at  Re  = 
100  leads  to  a  more  stable  flow  field  beyond  x/D  =  3. 
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Figure  4.9  Mean  flow  (a,b)  and  RMS  velocity  distributions  (c,d).  Left,  uncontrolled,  right, 
controlled  case.  The  cylinder  is  centered  at  (0,0)  and  of  diameter  1,  flow  from  left  to  right. 
Negative  isocontours  are  dashed,  positive  isocontours  are  solid  lines. 

Discussion 

f 

We  used  Proper  Orthogonal  Decomposition  (POD)  as  a  tool  to  process  multiple  sensor 
signals  into  a  global  estimate  of  the  flow  state.  POD  allows  for  stable  global  wake  state 
estimation,  enables  multi  sensor  evaluation  and  eliminates  artifacts  of  local  sensing,  i.e.  sensing 
at  nodes  of  the  vortex  street.  It  also  allows  for  an  accurate  state  estimate  when  the  effect  of  the 
controller  causes  major  changes  both  in  the  mean  flow  and  the  rms  amplitudes  of  the  fluctuating 
velocity  components.  However,  we  find  it  necessary  to  account  for  the  changes  in  the  mean  flow 
by  adding  a  mean  flow  mode  (mode  5)  to  the  model. 
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Figure  4.10  Linear  stability  analysis  of  unforced  and  feedback  controlled  flow  fields 


While  we  used  only  Mode  1  for  closing  the  feedback  loop,  all  the  higher  order  POD 
modes  experienced  proportional  reductions  in  mode  amplitude.  This  suggests  a  strong  coupling 
between  all  modes,  and  implies  that  the  existence  of  the  higher  order  modes  is  conditional  on  the 
presence  of  the  fundamental  modes.  This  confirms  the  results  of  our  previous  work. 

While  feedback  control  was  able  to  stabilize  the  near  wake  of  the  cylinder,  vortex 
formation  still  occurred,  though  further  downstream.  While  the  reasons  for  this  are  not  entirely 
clear,  we  suggest  several  possible  causes.  The  change  in  the  mean  flow  caused  by  the  controller 
lengthens  the  recirculation  zone.  This  moves  the  vortex  formation  location  further  downstream 
and  causes  a  reduction  in  both  drag  and  rms  lift  force.  While  both  of  these  effects  are  desired,  the 
downstream  shift  in  vortex  formation  location  causes  a  larger  spatial  separation  between  the 
actuation,  which  remains  at  the  cylinder,  and  the  oscillations  the  actuator  attempts  to  cancel.  This 
requires  both  more  actuation  input,  and  also  an  adjustment  in  the  actuation  phase  in  order  to 
account  for  the  time  a  given  disturbance  takes  to  travel  from  the  actuator  to  the  vertex  formation 
location.  At  the  same  time  the  disturbances  caused  by  the  actuator  need  to  travel  through  a  region 
of  the  flow  which,  while  stabilized,  is  only  stabilized  within  a  narrow  range  of  phase  angles.  If 
the  far  wake  requires  a  phase  angle  for  stabilization  that  at  the  same  time  destabilizes  the  near 
wake,  a  physical  limit  has  been  reached  in  terms  of  what  can  be  achieved  given  the  actuator 
location.  This  effect  may  limit  the  spatial  range  for  which  stabilization  can  be  achieved  with  the 
current  actuator  setup.  This  is  particularly  true  since  the  vortex  suppression  achieved  in  this 
investigation  decreases  the  stability  of  the  feedback  controlled  flow  field  and  extends  the 
unstable  region  of  the  flow  further  downstream  of  the  cylinder,  as  was  demonstrated  by  the  linear 
stability  investigations  shown  in  the  previous  section.  ,  <■ 

Despite  all  these  problems,  we  were  able  to  suppress  the  oscillations  in  the  near  wake 
without  actively  modifying  the  mean  flow  or  changing  the  separation  point  using,  for  example, 
momentum  injection.  Thus  this  effort  shows  that  the  cylinder  wake  flow  can  be  improved  in 
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terms  of  drag  and  unsteady  lift  by  feedback  control.  For  this  reason,  one  would  expect  the  current 
control  approach  to  be  applicable  to  wake  flows  with  fixed  separation  points,  like  the  flow 
around  a  D-  shaped  cylinder.  The  same  cannot  be  said  for  approaches  that  aim  at  moving  the 
separation  point  aft,  for  example  by  tripping  the  boundary  layer  or  using  blowing  and  suction 
upstream  of  the  separation  point  to  delay  separation. 

In  parallel  with  this  paper,  Gerhard  et  al.  also  employ  a  low  order  model  based  approach 
to  control  the  circular  cylinder  wake.  While  details  of  the  low  order  model,  estimation  technique 
and  control  algorithm  employed  are  different  from  our  approach,  the  qualitative  effects  of  the 
control  on  the  flow  field  appear  to  be  very  similar.  This  indicates  that  the  limitations  encountered 
in  both  of  our  control  approaches  may  be  inherent  to  actuation  authority  -  both  approaches  use  a 
localized  actuation  in  a  small  portion  of  the  flow  field  -  and  inherent  to  the  physics  of  the  flow 
field  itself. 

Overall,  we  were  able  to  reduce  the  effect  of  vortex  shedding  on  both  the  unsteady  lift 
and  the  vortex  induced  portion  of  the  pressure  drag  by  about  an  order  of  magnitude. 

4.6  Circular  Cylinder:  3-D  High  Reynolds  Number  Baseline  CFD 

Circular  Cylinder  calculations  at  a  subcritical  Reynolds  number  of  3900  have  been 
conducted  by  Hansen  and  Forsythe  (2003).  Although  this  work  was  not  performed  as  part  of  this 
STTR,  it  is  included  (briefly)  in  this  report  since  it  forms  the  baseline  for  similar  calculations  of 
the  “D”-shaped  cylinder  and  NACA  0015.  A  comprehensive  set  of  experiments  and  other 
computations  are  available  at  this  Reynolds  number  and  were  used  for  validation  of  the  solution 
method.  At  this  Reynolds  number,  the  boundary  layers  are  laminar,  and  the  shear  layers 
transition  to  turbulent  in  the  wake.  For  the  computations,  Large  Eddy  Simulation  (LES)  was 
applied  with  no  explicit  subgrid  scale  model.  This  means  that  the  numerical  dissipation  of  the 
solver  was  used  to  model  the  subgrid  scale  stresses.  LES  was  used  rather  than  DES  since  at  this 
Reynolds  number  the  boundary  layers  are  laminar,  so  a  hybrid  RANS/LES  method  such  as  DES 
holds  no  advantage  over  LES.  LES  assumes  sufficient  grid  resolution  to  resolve  the  large  scale 
energy  containing  eddies.  A  large  part  of  the  study  was  a  comprehensive  grid  refinement  study 
to  determine  the  necessary  grid  resolution.  Starting  from  a  baseline  grid,  two  other  grids  wree 
created.  Each  grid  refinement  was  achieved  by  increasing  the  cell  density  by  a  factor  of  V2  in 
each  coordinate  direction,  resulting  in  volume  cell  counts  of  442,018,  1,230, 7 10, -and  3,258,000, 
for  grids  A,  B,  and  C,  respectively.  The  effect  of  grid  refinement  on  the  flow  is  shown  in  Figure 
4. 1 1,  with  increasing  grid  resolution  resolving  more  structure. 
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Figure  4.11.  Top  and  side  views  of  isosurfaces  of  non-dimensional  spanwise  velocity 
perturbations  v'/U„  =  ±0.2 .  Red  shows  the  positive  valued  isosurfaces.  Grid  A  results  are  on  top, 
grid  B  results  are  below. 

Figure  4.12  shows  the  coefficient  of  pressure  on  the  cylinder  surface  for  the 
different  grids  used.  The  experimental  data  of  Norberg  (1987)  at  Re=3000  and  Re=4020  provide 
pressure  data  over  the  entire  surface.  Only  the  coefficient  of  backpressure  (Cp  at  0=180)  is 
available  for  Re=3900.  The  grid  C  results  match  very  closely  with  the  Re=3000  experiment, 
particularly  beyond  78  degrees.  The  actual  Re=3900  pressure  coefficient  line  is  judged  to  lie 
somewhere  close  to  midway  between  the  grid  B  and  grid  C  results  based  on  the  Re=3900  base 
pressure  coefficient.  Both  grids  B  and  C  seem  to  be  of  sufficient  resolution,  while  grid  A  is 
clearly  under-resolved. 

Hansen  and  Forsythe  (2003)  also  examined  wake  statistics  such  as  mean  centerline 
velocity  and  Reynolds  stresses,  and  for  the  two  finer  grid  there  was  good  agreement  with  the 
experiments. 


44 

STTR  Final  Report,  July  2004 


Cobalt  Solution,  LLC 
USAF  Academy 
STTR  Topic  AF03T007 
Contract:  F49620-03-C-0108 


Figure  4.12.  Cylinder  Surface  Coefficient  of  Pressure.  Experiment  is  from  Norberg . 

grid  A; - grid  B;  —  grid  C;  X  Tremblay  (2000)  DNS;  ■  experiment,  Re=  3000; 

▼experiment,  Re=3900;  O  experiment,  Re  =  4020. 


4.7  Summary 

The  effect  of  feedback  flow  control  on  the  wake  of  a  circular  cylinder  at  a  Reynolds 
number  of  100  is  investigated  in  direct  numerical  simulation.  Our  control  approach  uses  a  low 
dimensional  model  based  on  proper  orthogonal  decomposition  (POD).  The  controller  applies 
linear  proportional  and  differential  feedback  to  the  estimate  of  the  first  POD  mode.  The  range  of 
validity  of  the  POD  model  is  explored  in  detail.  Actuation  is  implemented  as  displacement  of  the 
cylinder  normal  to  the  flow.  We  demonstrate  that  the  threshold  peak  amplitude  below  which  Ijie 
control  actuation  ceases  to  be  effective  is  in  the  order  of  5%  of  the  cylinder  diameter. 

The  closed  loop  feedback  simulations  explore  the  effect  of  both  fixed  phase  and  variable 
phase  feedback  on  the  wake.  While  fixed  phase  feedback  is  effective  in  redrfcing  drag  and 
unsteady  lift,  it  fails  to  stabilize  this  state  once  the  low  drag  state  has  been  reached.  Variable 
phase  feedback,  however,  achieves  the  same  drag  and  unsteady  lift  reductions  while  being  able 
to  stabilize  the  flow  in  the  low  drag  state.  In  the  low  drag  state,  the  near  wake  is  entirely  steady, 
while  the  far  wake  exhibits  vortex  shedding  at  a  reduced  intensity.  We  achieved  a  drag  reduction 
of  15%  of  the  drag,  and  lowered  the  unsteady  lift  force  by  90%.  3-D  calculations  at  a  higher 
Reynolds  number  have  been  performed  and  validated,  and  are  available  for  subsequent  control 
attempts  in  phase  II. 
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5.  Application  II:  “D”  Shaped  Cylinder  Wake  Control 

5.1  “D”  Shaped  Cylinder:  Application  Background 

In  order  to  demonstrate  the  applicability  of  the  toolbox  to  new  flow  geometries,  the  D 
shaped  cylinder  geometry  was  considered.  The  main  motivation  to  include  the  D-shaped  cylinder 
in  the  AFC  project  is  because  it  has  a  fixed  separation  point.  This  means  that  the  only  way  to 
suppress  the  vortex  shedding  is  by  means  of  controlling  the  absolute  instability  of  the  wake  and 
not  by  moving  the  separation  point.  Additionally,  due  to  the  fixed  separation  point  this  geometry 
is  an  excellent  test  bed  for  demonstrating  actuation  using  a  blowing  and  suction  slot. 

CFD  computations  have  been  performed  on  a  grid  developed  for  the  “D”  shaped  cylinder. 
While  we  initially  used  a  Reynolds  number  of  150  to  make  sure  that  the  flow  is  two-dimensional 
and  laminar,  we  found  that  due  to  the  thick  boundary  layers  the  D  shaped  cylinder  flow  was 
much  more  stable  than  the  circular  cylinder  flow  at  the  same  Reynolds  number.  This  extends  the 
computation  times  for  both  open  loop  forced  and  feedback  controlled  simulations.  A  flow 
visualization  study  in  the  water  tunnel  revealed  that  the  critical  Reynolds  number  for  2D  to  3D 
transition  of  the  D  cylinder  flow  is  in  the  order  of  380. 

For  the  above  reasons,  it  was  decided  to  increase  the  Reynolds  number  to  300  for  the 
investigation  of  the  open  loop  forced  and  future  closed  loop  simulations.  Figure  5.0  compares  the 
two  different  Reynolds  numbers.  It  can  be  seen  that  the  higher  Reynolds  number  establishes  a 
limit  cycle  vortex  shedding  within  fewer  cycles.  This  speeds  up  the  computations,  since  the 
minimum  possible  time  step  is  a  fixed  fraction  of  a  shedding  cycle. 

•  Re=300 

•  M=0.1 

•  At*  =  0.005  shedding  cycles  (non-dimensional  time  step) 
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Figure  5.0  Lift  and  Drag  forces  during  the  startup  of  the  CFD  simulation.  Left,  Re  =  150;  right, 

Re  =  300 


5.2  “D”  Shaped  Cylinder:  Surface  Mounted  Sensor  Estimation 

Recent  research  on  closed-loop  control  of  the  von  K&rman  wake  instabilities  have 
addressed  the  issue  of  sensor  placement  and  number  based  on  non-intrusive  sensors  in  the  wake. 
This  approach  may  not  always  be  implemented  and  it  is  important  to  develop  an  effective 
method  for  sensor  placement  and  number  based  on  body  mounted  sensors.  These  sensors  may 
measure  skin  friction  or  surface  pressures,  as  done  in  this  effort.  The  advantages  of  surface 
mounted  sensors  are: 

•  Simple,  relatively  inexpensive  and  reliable 

•  Essential  for  real-life,  closed-loop  flow  control  applications  where  the  direct 
measurement  of  the  wake  flow  field  is  cumbersome  (if  not  impossible) 

•  Enable  “nearly  collocated”  sensors  and  actuators,  which  eliminates  substantial  phase 
changes  (affects  controller  design) 

In  this  section,  we  demonstrate  and  validate  a  systematic  approach,  based  on  the  spatial 
Eigen-functions  of  the  POD  model,  for  determining  sensor  number  and  location  for  estimation  of 
the  truncated  POD  states  of  a  cylinder  wake.  A  simple  model  of  the  flow-field  was  sought  to 
design  sensor  configuration  for  feedback  control  algorithms.  Numerical  simulations  were 
conducted  on  COBALT  solver  V.2.02.  In  this  effort,  the  solver  was  used  for  direct  numerical 
solution  of  the  Navier  Stokes  equations  with  second  order  accuracy  in  time  and  space.  An 
unstructured  two-dimensional  grid  with  120,000  nodes  and  115,000  elements  was  used.  The 
cylinder  geometry  comprises  of  a  semi  ellipse  with  base  height  of  7n1lm  and  length  of  61.25mm. 
Additional  simulation  parameters  are  as  follows: 

•  Initial  perturbation:  AOA=l  °  (to  kick  off  the 
vortex  shedding) 
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•  Damping  Coefficients: 

o  Advection  =  0.01 
o  Diffusion=  0.00 

•  Reynolds  Number:  Re  =  300 

•  Free  stream  velocity:  Ujnr=  34  [m/s] 

•  Shedding  Frequency  8.02  Hz 

=>  St  =  0.165 

•  Time-step  =  0.5  ms 

=>  250  time  steps  per  shedding  cycle 

The  simulation  was  triggered  by  skewing  the  incoming  mean  flow  by  a  =  1  degrees  to 
introduce  an  initial  perturbation.  The  validity  of  the  CFD  model,  used  in  this  effort  was 
established  comparing  the  Strouhal  number  from  the  simulation  to  a  water  tunnel  experiment. 
The  water  tunnel  experiment  yielded  St  =  0.163,  therefore  there  is  good  agreement  since  the 
difference  between  simulation  and  experiment  is  far  less  than  the  experimental  uncertainty  of 
0.01. 

In  this  effort,  the  flow  field  in  the  wake,  represented  by  the  Pressure  and  stream-wise 
Velocity  variables  are  obtained  from  the  numerical  solution  of  the  Navier-Stokes  equations 
obtained  using  the  CFD  model.  All  the  100  snap-shots  (~  15  shedding  cycles)  were  equally 
spaced  in  time.  The  snap-shots  were  taken  after  ensuring  that  the  cylinder  wake  reached  steady 
state.  For  control  design  purposes,  the  POD  method  enables  the  Navier-Stokes  equations  to  be 
modeled  as  a  set  of  ordinary  differential  equations  (O.D.E.).  At  first,  the  flow  field  data  is  loaded 
and  arranged  from  the  CFD  solver  of  the  "D"  shaped  cylinder  wake  at  Re  =  300.  The 
decomposition  of  this  component  of  the  velocity  field  (a  similar  representation  may  be  made  for 
the  Pressure  field)  is  as  follows: 

u(x,y,t)  =  U(x,y)  +  u(x,y,t)  (5.1) 


where  U  denotes  the  mean  flow  and  u  is  the  fluctuating  component  that  may  be  expanded  as: 


u(x,y,t)  =  2 ak  (t)0lk\x,y)  (5.2) 

*= i  ' 

where  a*(t)  denotes  the  time-dependent  coefficients  and  <pi(x,y)  represents  the  non- 
dimensional  spatial  Eigen-functions  of  the  Velocity  (see  Fig.  5.1)  and  Pressure  (see  Fig.  5.2) 
determined  from  the  POD  procedure.  From  an  ensemble  of  snapshots,  the  'mean  snapshot'  is 
computed  and  then  this  mean  is  subtracted  from  each  member  of  the  ensemble/  This  is  done 
primarily  for  reasons  of  scale;  i.e.  the  deviations  from  the  mean  contain  information  of  interest 
but  may  be  small  compared  with  the  original  signal. 

Next,  the  empirical  correlation  matrix  is  computed  using  the  inner  product.  Solving  the 
Eigen-value  problem,  the  Eigen-values  and  the  orthogonal  spatial  Eigen-functions,  <pi(x,y),  are 
obtained.  Since  the  Eigen- values  measure  the  relative  energy  of  the  system  dynamics  contained 
in  that  particular  mode,  they  may  be  normalized  to  correspond  to  a  percentage  of  the  system 
energy.  For  the  current  working  point  (Re=300),  the  Eigen-values  for  the  Pressure  and  U- 

48 

STTR  Final  Report,  July  2004 


lotions 


Cobalt  Solution,  LLC 
USAF  Academy 
STTR  Topic  AF03T007 
Contract:  F49620-03-C-0108 


Velocity  are  presented  in  Table  5.1.  Note  that  the  great  majority  of  the  energy,  associated  with 
the  POD  procedure,  is  located  in  the  first  two  modes. 


Pressure 

Velocity  (U) 

Mode 

POD 

POD 

Eigen-Values 

Eigen-Values 

[%] 

I%1 

H 

1 

49.31 

43.14 

Eft 

4.73 

2.92 

El 

4.68 

2.85 

Total 

98.92 

98.22 

Table  5.1:  Eigen-Values  of  Wake  Flow-Field 


Finally,  the  time  histories  of  the  temporal  coefficients  of  the  POD  model,  cik(t),  are 
determined  using  the  extracted  spatial  modes  and  the  snapshots  of  the  unforced  flow.  However,  it 
is  not  possible  to  obtain  a  direct  measurement  of  a*,  which  is  why  it  needs  to  be  estimated  from 
direct  measurements  such  as  body  mounted  sensors.  An  important  aspect  of  reduced  order 
modeling  concerns  truncation.  How  many  modes  are  important  and  what  are  the  criteria  for 
effective  truncation? 


Model  Mode2 

2 
0 

-2 

0  5  10  0  5  10 


Mode3  Mode4 


Figure  5.1.  Eigen-Functions  of  stream-wise  Velocity  of  the  “D”-shaped  cylinder  Wake  Flow 
Field  Solid  lines  are  positive,  dashed  lines  are  negative  isocontours 
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Figure  5.2  Eigen-Functions  of  Pressure  of  wake  Flow  Field 
Solid  lines  are  positive,  dashed  lines  are  negative  isocontours 


The  answers  to  the  above  questions  have  been  addressed  by  Cohen  et  al  (2003).  That 
effort  demonstrated  that  control  of  the  POD  model  of  the  von  Karman  vortex  street  in  the  wake 
of  a  circular  cylinder  at  low  Reynolds  numbers  (Re~100)  is  enabled  using  just  the  first  mode. 
Furthermore,  feedback  based  on  the  first  mode  alone  significantly  attenuated  all  the  other  modes 
in  the  four-mode  POD  model. 

In  view  of  the  above  result,  in  this  effort,  truncation  of  the  POD  model  takes  place  after 
the  first  three  modes,  which  contain  95-97%  of  the  energy,  as  seen  in  Table  5.1.  At  this  point,  it 
is  imperative  to  note  the  difference  between  the  number  of  modes  required  to  reconstruct  the 
flow  and  the  number  of  modes  required  for  effective  low-dimensional  modeling  for  control 
design  purpose.  In  this  effort,  we  are  interested  in  estimating  only  those  modes  required  for 
closed-loop  flow  control.  On  the  other  hand,  a  more  accurate  reconstruction  of  the  velocity  and 
pressure  field,  based  on  a  low-dimensional  model,  may  be  obtained  using  betweerf  4-8  modes. 

The  quintessential  question  is  whether  an  effective  estimate  of  the  states,  of  the  3  mode 
low-dimensional  model  coefficients,  a*,  can  be  estimated  based  on  body  mounted  pressure 
sensors.  The  answer  is  positive  and  the  details  that  provides  the  estimate  of  the  first  three  modes, 
a /-03,  are  presented  in  the  next  section. 

The  time  histories  of  the  temporal  coefficients  of  the  flow  field  Velocity  and  Pressure 
variables  are  determined  by  introducing  the  POD  spatial  Eigen-functions  into  the  flow  field  data 
using  the  least  squares  technique.  The  intent  of  the  proposed  strategy  is  that  the  pressure 
measurements  provided  by  the  body  mounted  pressure  sensots  are  processed  by  the  estimator  to 
provide  the  estimates  of  the  first  three  temporal  modes  of  the  flow  field  variables  stream-wise 
Velocity  and  Pressure.  The  estimation  scheme,  based  on  the  linear  stochastic  estimation  (LSE) 
procedure  introduced  by  Adrian  (1977),  predicts  the  temporal  amplitudes  of  the  first  three  POD 
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modes  from  a  finite  set  of  pressure  measurements  obtained  from  the  CFD  solution  of  the 
uncontrolled  cylinder  wake.  All  the  measurements  were  taken  after  ensuring  that  the  cylinder 
wake  flow  regime  converges  to  steady  state  vortex  shedding.  The  mode  amplitudes,  ay-aj,  will 
be  mapped  onto  the  extracted  sensor  signals  from  the  pressure  sensors,  us,  as  follows: 

«.(o=i;c>,(')  <5-3> 

.9=1 

where  m  is  the  number  of  sensors  and  Cs  represents  the  coefficients  of  the  linear  mapping.  The 
effectiveness  of  a  linear  mapping  between  for  velocity  measurements  and  POD  states  has  been 
experimentally  validated  by  Cameron  et  al  (2004).  The  coefficients  C”,  (n  =1-3;  s  =  1,  m)  in 
Equation  (5^)  are  obtained  via  the  linear  stochastic  estimation  method  from  the  set  of  discrete 
sensor  signals  and  temporal  mode  coefficients,  aj  -  03. 

The  issue  of  sensor  placement  and  number  has  been  dealt  with  in  an  ad-hoc  manner  in 
published  studies  concerning  closed-loop  flow  control.  For  effective  closed-loop  control  system, 
the  following  questions  need  to  be  answered: 

•  How  many  sensors  are  required? 

•  Where  should  the  sensors  be  placed? 

•  What  are  the  criteria  forjudging  an  effective  sensor  configuration? 

•  What  are  the  robustness  characteristics  of  a  given  sensor  configuration? 

In  this  effort,  an  attempt  will  be  made  to  emulate  some  of  the  proven  successes  from  the 
field  of  structural  control.  Heuristically  speaking,  when  some  very  fine  dust  particles  are  placed 
on  a  flexible  plate,  excited  at  one  of  its  natural  frequencies,  after  a  short  while  the  particles 
arrange  themselves  in  a  certain  pattern  typical  of  those  frequencies.  The  particles  will  be 
concentrated  in  the  areas  that  do  not  experience  any  motion  (the  nodes).  On  the  other  hand,  at  the 
areas  where  the  motion  is  large  (the  intemodes)  will  be  clean  of  particles.  It  is  at  the  intemodes 
that  the  vibrational  energy  of  that  particular  mode  is  at  a  maximum  and  sensors  placed  at  these 
locations  are  extremely  effective  in  estimating  that  particular  mode. 

The  above  heuristic  approach  has  been  used  by  Siegel  at  al  (2003)  in  locating  effective 
sensor  placement  for  acceleration  feedback  control  to  alleviate  tail  buffeting  of  a  high 
performance  twin  tail  aircraft.  Note  the  usage  of  the  term  “effective  sensor  configuration”  as  it  is 
based  on  validated  heuristicsy  as  opposed  to  “optimal  sensor  configuration”  that  results  from  a 
mathematically  optimal  pattern  search  for  a  sensor  configuration.  So,  what  needs  to  be  done  to 
determine  an  effective  sensor  configuration  is  to  find  the  areas  of  energetic  modal  activity. 

CFD  data  provided  simulated  pressure  signals  at  286  locations  distributed  at 
equidistances  all  along  the  surface  of  the  cylinder  (see  the  cylinder  surface  grid  in  Fig.  5.3).  A 
three-  step  procedure  is  proposed  for  determining  sensor  placement  and  number  as  follows: 

•  Run  the  POD  procedure  on  the  286  pressure  signals.  Examine  the  frequencies  of  the 
energetic  modes.  In  this  effort,  three  energetic  modes  are  found  and  their  frequencies 
correspond  to  the  fundamental  von  Karman  shedding  frequency  (first  two  modes)  and 
the  next  higher  frequency  (third  mode). 
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The  spatial  Eigen-functions  obtained  from  the  POD  procedure  provides  information 
concerning  the  locations  where  the  modal  activity  is  at  its  highest  (see  Fig.  5.5). 
Examine  the  maxima/minima  of  the  spatial  Eigen-functions  of  the  286  Pressure 
signals. 

Place  sensors  at  the  energetic  maxima  and  minima  of  each  mode  as  shown  in  Fig.  5.3. 
The  locations  of  the  sensors  in  Table  5.2  are  referenced  in  terms  of  the  coordinates, 
non-dimensionalized  with  the  model  base  height  H,  namely,  X/H  and  Y/H.  Three 
sensors  are  placed  at  the  upper  surface  of  the  cylinder  near  the  trailing  edge  and  three 
more  sensors  placed  symmetrically  at  the  lower  surface.  These  six  sensors  target 
modes  1  and  2.  Finally,  four  sensors  are  located  at  the  base  of  the  cylinder,  targeting 
mode  3.  The  time  histories  of  the  10  sensors  are  presented  in  Fig.  5.4.  Note  the  two 
distim^  frequencies  picked  up  by  the  sensors. 

Run  the  LSE  procedure  to  obtain  the  transformation  matrix  Cs  and  to  obtain  the 
estimates  of  the  U-Velocity/Pressure  POD  mode  coefficients.  The  estimated  versus 
desired  mode  amplitude  plot,  for  the  above  sensor  configuration  is  presented  in  Fig. 
5.6.  The  estimates  resulting  from  the  10  sensor  configuration  are  very  accurate  as  seen 
in  Table  5.2. 


For  convenience,  this  RMS  error  (defined  as  the  difference  between  the  RMS  of  the 
desired  extracted  mode  amplitudes  and  the  estimates  obtained  form  the  LSE  procedure)  is 
normalized  with  the  RMS  of  the  desired  extracted  mode  amplitudes,  presented  as  a  percentage. 
The  resulting  error  and  the  number  of  sensors  may  be  integrated  together  into  a  cost  function. 
The  purpose  of  the  design  process  would  then  be  to  select  the  configuration  that  minimizes  this 
cost  function.  For  this  configuration,  the  RMS  estimation  errors  are  provided  in  Table  5.2. 
Considering  the  fact  that  the  LSE  is  providing  U-Velocity/Pressure  estimates  of  the  wake  flow 
field,  the  RMS  values  in  Table  5.2  are  low  and  can  be  used  for  closed-loop  flow  control  using  a 
moderately  robust  controller.  Also,  provided  in  Table  5.2  are  the  RMS  errors  for  sensor 
configurations  consisting  of  2,  4  and  6  sensors.  Note  that  the  three  latter  configurations  enable 
the  estimation  of  modes  1  and  2  quite  well  but  not  of  the  3rd  since  they  do  not  include  a  sensor 
that  targets  the  higher  frequency. 


t 


Number 

of 

Sensors 

Sensor  Locations 

(x/H,y/H) 

Mode  1 

RMS  Error 
[%] 

Mode  2 

RMS  Error 

„  i%i 

Mode  3 

RMS  Error  [%] 

POD  Coeff. 

Estimated 

,  (in  Wake) 

2 

(-0.07,  0.50),  (-0.10,-0.50) 

22.74 

60.27 

Only  Modes  1  &  2 
Targeted 

Pressure 

■ 

2.99 

8.07 

Only  Modes  1  &  2 
Targeted 

Pressure 

6 

(-0.07,0.50),  (-0.10,-0.50) 
(-0.10,0.50),  (-0.07,-0.50) 
(-0.13,0.50),  (-0.13,-0.50) 

0.78 

2.17 

Only  Modes  1  &  2 
Targeted 

Pressure 
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10 

(-0.07, 0.50),  (-0.10,-0.50) 
(-0.10, 0.50),  (-0.07,-0.50) 
(-0.13,  0.50),  (-0.13,-0.50) 
(0,  0.42),  (0,-0.01) 

(0,  0.24),  (0,0.08) 

0.77 

2.10 

20.90 

Pressure 

10 

Same  10  sensor 
configuration  as  above 

2.30 

2.30 

12.5 

UJVelocity 

Table  5.2:  Sensor  Coordinates  and  RMS  Estimation 
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Time  Signals  of  Pressure  Taps 


18.52 


2.5 


2.6  2.7  2.8  2.9  3  3.1  3.2  3.3  3.4 

Time  [s] 

Figure  5.4.  Time  Histories  of  Pressure  taps  for  the  10  Sensors 
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Eigen-Functions  along  X  axis 


Figure  5.5.  Spatial  Surface  Eigen-functions  for  286  Pressure  sensors  (note  that  the  mode 
amplitudes  are  plotted  against  x/H  -  in  the  top  Figure;  and  against  Y/H  in  the  bottom  Figure) 
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Comparison  Mode  Amplitudes  based  on  Sensors  (*)  and  Full  Flow  Field  (-) 


Comparison  Mode  Amplitudes  based  on  Sensors  (*)  and  Full  Flow  Field  (-) 


3.5 


Figure  5.6.  Estimation  of  POD  Time  Coefficients  for  10  Sensor  configuration 
(UJVelocity  -  Top  Figure;  Pressure  -  Bottom  Figure) 
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5.3  “D”  Shaped  Cylinder:  Blowing/Suction  Modeling 

Figure  5.7  shows  the  grid  implementation  of  blowing  and  suction  slots  into  the  D 
shaped  cylinder  geometry.  The  actuators  are  modeled  by  short  channels  at  an  angle  of  30 
degrees  to  the  mean  flow.  At  the  end  of  the  channel  a  either  periodic  of  feedback 
controlled  velocity  is  prescribed,  leading  to  an  unsteady  flow  through  the  channel  into  the 
flow  field. 


Figure  5.7  Computational  grids  developed  for  modified  D-shaped  cylinder 
using  blowing  and  suction  actuation. 

Figure  5.8  demonstrates  the  effect  of  periodic  open  loop  blowing  and  suction 
through  these  patches  at  the  bottom  of  short  ducts  at  an  angle  of  30  degrees  to  the 
ffeestream  direction.  The  forcing  frequency  is  equal  to  the  natural  shedding  frequency, 
while  the  peak  velocity  during  a  sinusoidal  blowing  and  suction  cycle  is  12.5%  or  100% 
of  the  ffeestream  velocity,  respectively.  The  top  and  bottom  slot  were  operated  180 
degrees  out  of  phase  for  all  investigations.  The  plots  show  isocontours  df  vorticity, 
demonstrating  the  local  effect  of  the  forcing.  The  blowing  and  suction  velocity  is  not 
large  enough  for  the  resulting  jet  to  penetrate  the  boundary  layer  in  either  case.  However, 
with  the  higher  velocity  more  vorticity  is  injected  into  the  boundary  layer,  effectively 
cutting  of  the  forming  vortices.  Additionally,  the  vorticity  ejected  from  the  slots  then 
merges  with  the  vortices  of  the  Karman  vortex  street. 
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Figure  5.8  Instantaneous  snapshots  of  vorticity  during  a  blowing  and  suction  cycle.  Left, 
12.5%  of  the  freestream  velocity,  right,  100%  of  the  freestream  velocity. 

As  can  be  seen  in  the  transient  lift  and  drag  force  plots  shown  in  Figure  5.9,  the 
higher  velocity  forcing  creates  a  transient  behavior  that  is  quite  different  in  that  there  is 
significant  overshoot  before  the  wake  settles  into  a  phase  locked  state.  This  is  a  result  of 
the  vorticity  ejected  from  the  blowing  and  suction  slot  merging  with  the  vortex  street,  an 
effect  that  is  not  existent  for  smaller  velocities.  The  open  loop  investigations  used 
blowing  and  suction  peak  velocities  down  to  3%  of  the  freestream  velocity,  however,  for 
values  smaller  than  about  10%  the  transient  settling  times  increased  strongly.  This 
indicates  a  loss  of  controllability,  and  10%  was  determined  to  be  the  minimum  desirable 
actuation  amplitude  for  the  subsequent  closed  loop  controlled  simulations. 
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Figure  5.9  Transient  lift  (green)  and  drag  (blue)  forces  for  different  actuation  levels;  top 
left:  12.5%  of  the  freestream  velocity,  top  right:  25%,  bottom  left:  75%  and  bottom  right: 

100%  of  the  freestream  velocity. 

5  4  “D”  Shaped  Cylinder:  Feedback  Controlled  Results  and 
Discussion  , 

For  the  feedback  controlled  simulations,  the  spatial  modes  shown  in  figure  5.10 
were  employed  for  global  flow  state  estimation.  They  are  qualitatively  similar  to  the 
modes  found  in  the  circular  cylinder  wake  which  are  shown  in  Figure  4.3.  The  procedure 
used  to  develop  these  modes  is  identical  to  the  procedure  used  for  the  circular  cylinder, 
which  is  presented  in  the  previous  chapter.  The  sensor  arrangement  used  for  the  circular 
cylinder  consisting  of  a  rectangular  grid  of  35  sensors  was  also  kept. 
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Figure  5.10  POD  modes  used  for  feedback  state  estimation. 

A  scan  through  the  range  of  possible  controller  feedback  phases  showed  an 
increase  in  vortex  shedding  and  modal  amplitudes  for  phases  in  the  range  of  0  -  120 
degrees,  a  neutral  range  where  the  intensity  of  the  vortex  street  remained  about  the  same 
between  120  -  180  degrees,  and  reduction  in  vortex  shedding  strength  for  phase  angles 
beyond  180  degrees.  Figure  5.1 1  demonstrates  the  ability  to  reduce  vortex  shedding  with 
a  feedback  phase  of  210  degrees.  The  unsteady  lift  fluctuations  decrease  continuously 
over  about  20  shedding  cycles  and  are  stabilized  at  a  level  40%  lower  than  in  the  natural 
wake.  During  the  same  time  the  drag  force  decreases  by  about  10%.  This  behavior  is 
qualitatively  similar  to  the  circular  cylinder  wake,  while  the  reduction  levels  in  unsteady 
lift  and  mean  drag  force  have  not  yet  been  matched.  This  may  be  accomplished  by 
adjusting  the  feedback  parameters  further,  which  is  exceeded  the  scope  of  this  study. 
Figure  5.12  shows  the  changes  in  vortex  shedding  patterns  due  to  the  feedback.  The 
formation  length  has  increased,  and  the  vortices  being  shed  are  less  strong  than  in  the 
unforced  case.  These  findings  compare  well  with  the  circular  cylinder  results  as  well. 
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Figure  5.12  Unforced  (left)  and  feedback  controlled  (right)  vorticity  contours. 
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Figure  5.13  Unforced  (left)  and  feedback  controlled  (right)  mean  flow. 


The  transfer  of  energy  from  the  vortex  shedding  back  to  the  mean  flow  is  shown 
in  figure  5.13.  A  similar  lengthening  of  the  recirculation  zone  that  was  observed  in  the 
circular  cylinder  wake  is  achieved  in  the  D  shaped  cylinder  wake,  even  though  the 
increase  in  recirculation  zone  length  is  again  less  than  in  the  circular  cylinder  wake. 

While  the  geometry  of  the  circular  and  D  shaped  cylinders  are  very  different,  and 
different  actuation  methods  were  employed,  we  demonstrate  that  the  feedback  flow 
control  toolbox  is  well  suited  to  develop  and  test  control  algorithms  for  this  flow  field. 
The  approach  of  using  POD  for  global  flow  estimation  is  just  as  well  suited  for  the 
blowing  and  suction  slot  actuated  flow  field  as  it  is  for  the  translation  actuated  flow  field 
of  the  circular  cylinder.  Moreover,  the  response  of  both  flow  fields  to  the  feedback 
control  is  very  similar,  indicating  that  feedback  flow  control  is  applicable  to  different 
wake  flow  topologies.  Since  the  D  shaped  cylinder  features  fixed  separation  points,  the 
data  presented  shows  conclusively  that  our  control  approach  influences  the  vortex 
shedding  and  formation  in  the  wake  directly,  rather  than  modifying  the  stability  behavior 
of  the  mean  flow  by  using  separation  control. 

5.5  “D”  Shaped  Cylinder:  Experimental  Validation  and 
Discussion 

In  order  to  validate  the  computational  results,  the  vortex  shedding  frequency  of 
the  D  shaped  cylinder  simulations  were  compared  to  water  tunnel  experiments.  Figure 
5.14  shows  good  agreement  between  experiment  and  simulation,  indicating  sufficient 
resolution  of  the  numerical  model.  In  case  of  under  resolved  numerical  simulation  in 
unsteady  wake  flows  the  shedding  frequency  deviates  significantly,  as  has  been  shown 
for  circular  cylinder  wakes  in  literature. 
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Plot  of  Strouhal  Number  vs.  Reynolds  number  for  Water  tunnel  data  and  CFD  data 
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Figure  5.14  Strouhal  number  comparison  for  different  Reynolds  numbers. 
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5.6  “D”  Shaped  Cylinder:  3-D  High  Reynolds  Number  Baseline 
CFD 

Calculations  of  the  “D”  shaped  cylinder  were  performed  at  a  Reynolds  number 
based  on  base  height  of  3900,  using  the  experienced  gained  on  the  circular  cylinder.  The 
two-d  grid  used  for  the  closed  loop  control  cases  was  extruded  into  the  span  using  66 
points  and  4  base  diameters,  resulting  in  a  grid  with  7.5x1 06  cells.  This  was  larger  than 
the  circular  cylinder  grids  due  to  the  very  fine-spacing  in  the  x-y  plane.  LES  with  no 
explicit  model  was  applied  just  as  was  the  case  for  the  circular  cylinder.  A  non- 
dimensional  timestep  of  0.01  (made  non-dimensional  by  the  freestream  velocity  and  base 
height)  was  used.  A  visualization  of  the  flow  using  iso-surfaces  of  vorticity  is  shown  in 
Figure  5.15.  Note  that  the  shedding  is  occurring  at  different  phases  along  the  span, 
evidenced  by  spanwise  pressure  variations  near  the  back  end  of  the  cylinder.  This 
demonstrates  the  likely  necessity  of  spanwise  distributed  actuators.  The  wake  exhibits 
similar  structures  as  seen  for  the  circular  cylinder.  There  are  spanwise  rollers  connected 
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to  streamwise  counter-rotating  vortices.  These  results  will  be  available  for  use  in  phase 
II  for  3-D  attempts  at  control  of  a  higher  Reynolds  number  flow. 


Figure  5.15  LES  calculation  of  the  “D”  shaped  cylinder.  Isosurface  of  vorticity  colored 
by  pressure. 
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6.  Application  III:  NACA  0015  Airfoil  Separation  Control 

6.1  NACA  0015  Airfoil:  Application  Background 

Separation  control  on  the  NACA  0015  airfoil  serves  as  the  final  demonstration 
case  of  the  developed  flow  control  toolbox.  It  should  be  realized  that  in  this  case  the 
objective  of  flow  control  is  entirely  different  than  in  the  wake  cases  described  above. 
The  goal  here  is  the  delay  of  flow  separation,  and  it  is  clear  from  a  hydrodynamic 
stability  point  of  view  that  the  structures  in  the  boundary  layer  have  to  be  enhanced  to 
increase  the  momentum  transfer  from  the  free  stream  toward  the  wall.  This  is  essentially 
the  opposite  of  the  goal  in  control  of  the  wake,  where  flow  control  is  applied  to  reduce 
the  strength  of  the  vortices  in  the  wake.  However,  our  toolbox  was  successfully  applied 
to  this  case  as  well,  as  described  below. 

6.2  NACA  0015  Airfoil:  Computational  Model 

Calculations  were  performed  at  a  chord  Reynolds  number  of  Rec=l  0,000.  The 
grid  is  shown  in  figure  6.1.  As  a  first  step,  a  suitable  angle  of  attack  was  determined  by 
performing  unforced  computations  angles  of  attack  ranging  from  0°  to  1 5°.  From  these 
simulations,  it  was  determined  that  a=8°  is  optimal.  Once  a  test  case  was  defined,  we 
performed  open  loop  simulations  utilizing  the  blowing  and  suction  boundary  condition 
module  at  a  position  of  10%  chord  on  the  upper  airfoil  surface  to  verify  the  effectiveness 
of  blowing  and  suction  for  controlling  airfoil  separation.  The  simulations  confirmed  the 
applicability  of  periodic  forcing  for  separation  control.  Figure  6.2a  shows  a  snapshot  of 
the  instantaneous  pressure  field  of  the  unforced  case,  which  indicates  massive  separation 
near  the  leading  edge  and  the  development  of  large  vortical  structures  near  the  trailing 
edge  of  the  airfoil. 


Figure  6.1  Computational  grid  developed  for  NACA  0015 
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Figure  6.2:  Instantaneous  pressure  and  streamlines:  a)  unforced,  b)  forced 


6.3  NACA  0015  Airfoil:  Computational  Results  and  Discussion 

Scrutinizing  the  spatial  POD  modes  of  the  streamwise  velocity  component  of  the 
unforced  case  (Figure  6.3)  as  well  as  the  respective  temporal  coefficient  (Figure  6.4),  it 
becomes  clear  that  two  different  frequencies  are  present,  one  of  the  vortices  in  the 
separated  shear  layer  (Modes  1  and  2,  containing  96%  of  the  energy)  at  a  frequency  of 
F=25Hz  and  one  in  the  boundary  layer  (Modes  3  and  4,  containing  an  additional  3%)  at  a 
frequency  of  about  F=60FIz).  In  order  to  introduce  disturbances  in  the  most  efficient 
manner,  this  frequency  was  chosen  as  the  forcing  frequency. 


Model 


Mode2 


0.2  0.4  0.6  0.8 

Mode4 


Figure  6.3:  POD  modes  of  the  U  velocity  for  the  unforced  case 
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Figure  6.4:  Temporal  Mode  Coefficients  for  the  unforced  case 


Using  a  blowing  and  suction  slot  2.5%  chord  long  at  a  location  of  10%  chord  on 
the  upper  surface  of  the  airfoil,  periodic  blowing  and  suction  with  a  maximum  amplitude 
of  4m/s  (13%  of  the  ffeestream  velocity)  was  used  to  attempt  the  delay  of  separation. 
Figure  6.2b  shows  a  snapshot  of  the  instantaneous  flow  field  after  the  initial  transient. 
The  alternating  high  and  low  pressure  regions  indicate  the  presence  of  large  vortices  near 
the  airfoil  surface.  The  resulting  POD  modes  of  the  U  velocity  are  shown  in  Figure  6.1. 
It  can  be  seen  that  Modes  1  and  2  (containing  98%  of  the  energy)  of  the  forced  flow 
correspond  to  Modes  3  and  4  of  the  unforced  case.  Furthermore,  their  maxima  are 
significantly  closer  to  the  airfoil  surface,  indicating  a  reduction  of  size  of  the  separation. 
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Figure  6.1:  POD  modes  of  the  U  velocity.  Forced  case,  F=60Hz 

Concomitant  with  these  changes  in  flow  topology,  forcing  has  a  large  effect  on 
the  lift  and  drag.  Figure  6.  shows  the  time  history  of  the  lift  (Fy)  and  drag  (Fx)  for  the 
unforced  (left)  and  forced  (right)  cases.  Clearly,  active  control  not  only  markedly 
reduces  the  amplitude  of  the  force  fluctuations,  but  also  -  in  the  mean  -  increases  the  lift 
by  more  than  120%  and  reduced  the  drag  by  over  55%. 


t 
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a) 


b) 

Figure  6.5:  Lift  (left)  and  drag  (right)  for  a)  unforced  and  b)  forced  case 

t 

To  close  the  feedback  loop,  sensor  placement  has  been  studied  for  this  flow,  using 
the  heuristic  method  developed  earlier.  In  the  current  configuration,  20  sensors  are 
placed  at  the  relative  maxima  of  the  first  U  velocity  modes  for  both  the  forced  and 
unforced  cases  (see  Figure  6.6).  Using  the  time  signals  at  the  sensor  locations,  the  flow 
state  can  be  estimated  with  errors  of  less  than  10%.  Figure  6.7  shows  the  reconstruction 
of  the  time  signals  of  the  time  signals  of  the  unforced  flow  (left)  and  the  forced  flow 
(right). 
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**  Figure  6.6:  Grid  points  (blue)  and  sensor  location  (red) 


To  close  the  feedback  loop,  a  simple  proportional  feedback  controller  is  used  for 
the  forcing  amplitude  only.  In  order  to  determine  the  sensor  signal  best  suited  for 
feedback,  a  transient  simulation  was  performed.  In  this  simulation,  the  forcing  is  started 
during  the  calculation  to  observe  the  complete  transition  from  the  unforce'Id  (natural), 
separated  flow  to  the  forced.  Utilizing  the  POD  module  of  our  toolbox,  the  amplitude  of 
the  unforced  and  forced  mean  flow  modes  are  extracted  (see  figure  6.8).  As  can  be  seen 
in  the  figure,  the  flow  changes  from  the  unforced  state,  where  the  forced  mean  flow  mode 
amplitude  is  approx.  1  and  the  unforced  mode  amplitude  is  0,  to  the  forced  state  with  the 
amplitudes  reversed.  This  significant  change  in  mode  amplitudes  serves  as  a  good 
indication  of  the  flow  state,  and  with  an  appropriate  feedback  gain,  can  be  used  as  a 
feedback  signal.  However,  it  is  also  clear  that  a  minimum  threshold  amplitude  is  needed 
to  keep  the  flow  attached. 

In  the  simulation  presented  here,  the  amplitude  limits  were  set  to  Amax=5.5m/s  and 
Amin=2.5m/s.  With  these  limits,  lift  and  drag  as  a  function  of  time  are  shown  in  figure 
6.9.  When  compared  to  the  results  of  purely  periodic  blowing  and  suction  at  A=4m/s,  the 
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feedback  controlled  run  shows  slightly  lower  lift  and  higher  drag,  but  the  forcing 
amplitudes  are  considerably  lower.  Figure  6.10  shows  the  amplitude  of  the  unforced 
mean  flow  mode  as  a  function  of  time  (top)  and  the  accompanying  forcing  velocity 
(bottom).  Analyzing  figures  6.9  and  6.10,  it  becomes  clear  that  the  controller  is  able  to 
sense  the  deviation  of  the  sensor  signal  from  zero  and  accordingly  increase  the  forcing 
amplitude.  However,  because  of  the  time  delay  between  forcing  slot  and  the  sensors,  the 
controller  is  not  able  to  react  fast  enough  to  suppress  incipient  separation.  The  short 
bursts  of  separated  flow  lead  to  the  comparatively  large  temporal  fluctuations  of  lift  and 
drag. 

These  results  are  by  no  means  optimal,  but  clearly  demonstrates  that  our  toolbox 
can  be  applied  successfully  for  the  case  of  separation  control.  To  further  improve  the 
effectiveness  of  the  controller,  and  especially  to  take  the  time  delay  between  a  change  in 
actuation*  and  the  response  at  the  sensor  locations  into  account,  different  sensor  and 
controller  configurations  are  currently  under  investigation. 

>5f 


M«*n  flow  mode,  forced 


Mem  to*  mode,  unforced 
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Figure  6.8:  Temporal  development  of  mean  flow  modes.  Unforced  (blue)  and  forced  (red) 
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Figure  6.9:  Lift  (left)  and  drag  (right)  as  a  function  of  time  for  the  feedback  controlled  run 
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Figure  6.10:  Amplitude  of  unforced  mean  flow  mode  (top)  and  forcing  velocity  (bottom) 
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6.4  NACA  0015  Airfoil:  3-D  Baseline  CFD 


3-d  calculations  were  performed  for  the  NACA  0015  at  the  same  flow  conditions 
as  controlled  in  2-d,  i.e.  a  Reynolds  number  of  10,000,  and  an  angle  of  attack  of  8 
degrees.  The  2-d  grid  used  in  the  control  case  was  extruded  into  the  span  using  66  points 
and  one  chord  length,  leading  to  800,000  cells.  Simulations  were  performed  with  a  non- 
dimensional  timestep  of  0.01  (made  non-dimensional  by  the  chord  and  freestream 
velocity).  LES  with  no  explicit  subgrid  scale  model  was  applied  due  to  the  expected 
laminar  boundary  layers.  Flow  visualizations  are  show  in  Figure  6.1 1.  The  variation  in 
separation  in  the  spanwise  direction  is  seen,  again  emphasizing  the  need  to  progress  the 
closed  loop  tools  to  handle  3-D  cases.  The  resolved  content  in  this  simulation  is  fairly 
weak  compared  to  what  may  be  expected  at  this  Reynolds  number,  and  what  was  seen  for 
the  circular  and  D-shaped  cylinders.  Further  grid  refinement  may  be  required.  These 
results  will  be  available  for  subsequent  attempts  at  control  in  3-d  during  phase  II. 


Figure  6.1 1  LES  calculation  of  the  NACA  0015  airfoil.  Isosurface  of  vorticity  colored  by 
pressure.  Flow  shown  at  two  different  times. 
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7.  Summary  of  Phase  I  &  Recommendations  for  Phase  II 

At  the  end  of  Phase  I  of  this  project,  we  conclude  that  we  developed  a  toolbox 
that  meets  or  exceeds  the  capabilities  outlined  in  the  proposal.  We  were  able  to 
demonstrate  the  successful  application  of  the  toolbox  to  control  three  feedback  flow 
control  applications  as  proposed.  For  these  reasons,  in  comparing  our  recommendations 
for  phase  II  outlined  in  the  next  section,  you  will  find  them  to  be  identical  to  the  original 
suggestions  we  made  in  the  Phase  I  proposal,  with  more  detail  added. 

7.1  Main  Achievements  of  Phase  I 

The  mairifachievements  of  Phase  I  are  as  follows: 

* 

•  Toolbox  is  complete  and  fully  functional  and  tested 

•  Demonstrated  application  of  our  modeling  and  control  approach  to  three  different 
geometries,  namely  a  Circular  Cylinder,  D  shaped  Cylinder  and  NACA  0015 
airfoil 

•  Implemented  and  tested  feedback  controlled  boundary  conditions  in  the  Cobalt 
CFD  solver. 

•  Implemented  and  tested  hdf  output  in  the  Cobalt  CFD  solver  to  provide  data  to  the 
matlab  based  flow  control  toolbox  components 

•  Performed  3D  simulations  of  the  D-shaped  cylinder  and  the  NACA  0015  airfoil 

•  Overall,  we  showed  that  we  can  accomplish  low  dimensional  modeling,  control 
and  closed  loop  simulation  for  very  different  flow  fields  using  very  different 
actuation  approaches.  All  of  this  applies  to  2D  flow  fields. 

•  The  POD  algorithm  developed  deals  with  unstructured  data  and  delivers  accurate 
results  at  the  body  surface. 

•  A  shell  script  was  developed  that  starts  the  controller  automatically  with  the  CFD 
simulation  from  within  the  queue.  Thus  feedback  controlled  runs  can  be  queued 
just  like  any  other  CFD  simulation  run. 

These  achievements  meet  and  in  some  cases  exceed  the  goals  of  the  proposed 
research.  They  demonstrate  the  feasibility  of  the  approach  and  are  a  sound  basis  to  move 
ahead  into  Phase  II  of  the  STTR. 
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7.2  Recommended  Development  Approach  for  Phase  II 

1.  Extend  the  capabilities  of  the  Toolbox  to  3D  flow  fields 

Reasons  and  justification  to  do  this: 

•  Recently,  we  found  that  even  nominally  2D  flow  fields  like  the  cylinder  wake  turn 
3D  once  feedback  control  is  applied  (Siegel  et  al.  2004).  Furthermore,  higher  Re 
flows  tend  to  be  3D  anyways. 

•  This  suggests  that  it  in  order  to  provide  a  toolbox  that  enables  realistic  feedback 
flow  control,  it  is  imperative  to  expand  the  capabilities  to  the  3D  regime. 

•  V^Siile  some  parallel  research  teams  may  have  /  have  2D  closed  loop  capabilities, 
usually  they  use  home  brew  CFD  codes  that  cannot  easily  be  expanded  to  3D. 

•  Cobalt  has  excellent  3D  capabilities  with  DES  and  turbulence  modeling.  The 
proposed  plan  leverages  Cobalt’s  strength  in  LES  and  DES.  Cobalt’s  turbulence 
modeling  software  building  blocks  can  be  put  to  good  use  in  order  to  be  able  to 
build  low  dimensional  models  at  higher  Reynolds  numbers 

2.  Increase  the  Reynolds  number  into  realistic  Regimes 


Reasons  and  justification  to  do  this: 

•  While  we  have  learned  many  important  lessons  in  feedback  flow  control  working 
on  low  Reynolds  number  problems,  it  is  necessary  to  expand  the  Reynolds 
number  range  in  order  to  be  able  to  do  real  life  problems 

•  This  requires  the  3D  capabilities  outlined  above 

•  It  again  leverages  Cobalt’s  strength  in  LES  and  DES 

•  Cobalt’s  turbulence  modeling  software  building  blocks  can  be  put  to  good  use  in 
order  to  be  able  to  build  low  dimensional  models  at  higher  Re 

3.  Extend  the  capabilities  of  the  Toolbox  to  model  Body  Forces  for 
actuation 

Reasons  and  justification  to  do  this: 


•  Body  forces  are  a  generic  way  to  model  actuation.  They  can  be  implemented  at 
less  cost  in  a  CFD  simulation  than  using  a  high  resolution  grid  to  accurately 
resolve  actuator  geometry,  thus  improving  the  performance  of  the  simulation 

•  Body  forces  can  be  used  to  model  the  effect  of  plasma  actuators,  which  have 
shown  a  lot  of  promise  for  feedback  flow  control  applications  due  to  their  simple 
setup  without  any  moving  mechanical  parts 
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8.  Matlab  Modeling  and  Control  Toolbox  Suite 


This  chapter  is  intended  to  give  some  insight  into  the  functionality  of  the  main 
components  of  the  Matlab  based  feedback  flow  control  toolbox  developed  within  this 
program.  It  is  neither  complete,  nor  intended  to  be  a  manual  for  the  toolbox,  which  will 
be  developed  in  Phase  II  of  the  project.  It  lists  the  calling  convention  along  with  the 
documentation  of  the  modules,  which  gives  insight  into  both  functionality  and 
implementation  of  the  tools.  Note  that  none  of  the  many  sub-functions  are  listed,  but 
rather  the  functions  a  typical  user  will  interface  with.  Also  omitted  are  scripts  dedicated 
to  plotting  or  presentation  of  data.  The  overall  number  of  matlab  scripts  and  functions  in 
the  toolbox  specific  to  this  project  is  currently  in  the  order  of  250. 


function  [Energy-Content ,  U_Modes , varargout]  = 

POD_Modes (NModes ,  SubMean, varargin) 

% 

%  [EnergyContent ,  U_Modes]  =  POD_Modes (NModes ,  SubMean,  U_SnapShots) 

%  [EnergyContent,  U_Modes,  U_SnapShotMean]  =  POD_Modes (NModes ,  SubMean, 
U_SnapShots) 

%  [EnergyContent,  U_Modes,  V__Modes]  =  POD_Modes (NModes ,  SubMean, 
U_SnapShots,  V_SnapShots) 

%  [EnergyContent,  U_Modes,  U_SnapShotMean,  V_Modes,  V_SnapShotMean  3  = 
POD__Modes (NModes,  SubMean,  U_SnapShots,  V_SnapShots) 

% 

%  General  Inputs: 

%  NModes:  While  the  POD  yields  as  many  modes  as  there  are 

Snapshots,  only  NModes  will  be  returned  by  this  function 
%  SubMean:  SubMean  indicates  if  the  average  of  all  snapshots  is  to 

be  subtracted 

%  before  calculating  the  POD  modes: 

%  SubMean  =  1 :  Snapshot  Mean  WILL  be  subtracted 

%  SubMean  =  0:  Snapshot  Mean  WILL  NOT  be  subtracted 

% 

%  1  Dimensional  or  unstructured  Data: 
%+++++++++++++++++++++++++++++++++++++ 

%  The  u  and  (optional)  v  velocity  components  are  stored  in  matrices 
with  size  NxM .  t 

%  N  is  the  number  of  gridpoints  and  M  is  the  number  of  snapshots. 

%  The  Spatial  Modes  are  returned  in  U,V_Modes,  with  size  NxM. 

%  U, V_SnapShotMean  is  the  average  of  all  U,V  Snapshots,  dimension  is 

Nxl . 

% 

%  2  Dimensional  Data:  { 

%++++++++++++++++++++ 

%  The  2D  scalar  Data  Field (s)  is/are  in  a  3D  Matrix  called  U,V 
Snapshots,  size (Snapshots) = [Nx,  Ny,  M] 

%  where:  Nx  and  Ny  are  the  number  of  gridpoints  from  in  X  and  Y 
Direction  at  intervals  of  DeltaXY 

%  M  is  the  number  of  snapshots,  typically  at  different  instants  of  time 
%  The  spatial  Modes  are  returned  in  U,V  Modes,  with  sizes  [Nx  Ny 
NModes] 

% 

% 

%  The  Energy  Content  of  the  retained  modes  is  stored  in  EnergyContent. 
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% 

%  Revision  History:  Code  assembled  3-10-2004 


function  [TimeCoeff]  =  POD_Time_Coef f (U_SnapShots ,  ModesU, 
varargin) ; 

% [TimeCoeff]  =  POD_Time_Coef f (U_SnapShots ,  ModesU) 

% [TimeCoeff]  =  POD_Time_Coef f {U_SnapShots ,  ModesU,  MeanU) 

% [TimeCoeff]  =  POD_Time_Coef f (U_SnapShots ,  ModesU,  V_SnapShots,  ModesV) 

% [TimeCoeff]  =  POD_Time_Coef f (U_SnapShots ,  ModesU,  MeanU,  V_SnapShots, 
ModesV,  MeanV) 

%  Does  a  least  squares  fit  to  determine  the  temporal  coefficients  of  the 
POD  Modes 

%  U  Velocity  and  V  Velocity  as  well  as  ModesU  and  ModesV  need  to  have  the 
same  dimensions 

function  ReconData  =  POD_Reconstruct (Modes,  SnapShotMean, 

TimeCoeff) 

% 

%  ReconData  =  POD_Reconstruct (Modes ,  SnapShotMean,  TimeCoeff) 

%  Reconstructs  the  Snapshots  in  ReconData  from  Modes,  Mean  and  Time 
Coefficients 

%  5-4-04  SGS  supports  UNS  data 

function  OrthMatrix  =  POD_Orthogonality (Modes ,  DeltaXY) 

%  Checks  if  the  POD  Modes  are  orthogonal . 

%  If  they  are,  OrthMatrix  should  equal  the  Kronecker  delta  function 

function  [ObsMatrix,  EstimatedSignals,  MeanErrors, 
ErrorSignals]  = 

LinearStochasticEstimation (Observation,  ToBeEstimated) 


%************************************************************************ 
★  * 

%  Purpose  /  Usage 

%************************************************************************ 
★  * 

%  Linear  Mapping  of  Observations  onto  Estimated  States  or  vice  versa: 

% 

%  If  you  use  Time  signals  in  a  flow  field  as  Observations 
%  AND  POD  temporal  coefficients  as  ToBeEstimated  Signals  f 

%  you  will  develop  an  ObsMatrix  that  you  can  multiply  with  instantaneous 
%  Velocities  in  order  to  get  POD  Amplitude  estimates 

%  If  you  use  Time  signals  in  a  flow  field  as  ToBeEstimated  Signals 
%  AND  POD  temporal  coefficients  as  Observations 

%  you  will  develop  an  ObsMatrix  that~you  can  multiply  with  instantaneous 
%  Mode  Amplitudes  in  order  to  get  Velocity  Field  estimates 

%************************************************************************ 
*  * 

%  Data  Structures  /  Dimensions 

%**************************************************  *i*  ******************** 
** 

%  Observation  is  a  2D  matrix  size  M  x  N 
%  M  being  the  number  of  Observations 
%  N  being  the  number  of  time  steps  /  samples 
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%  ToBeEstimated  is  a  2D  matrix  that  contains  K  x  N  time  signals  measured 
in 

%  the  flow  field 

%  K  is  the  number  of  signals  to  be  estimated 
%  ObsMatrix  is  2D  and  has  K  x  M  Elements 

%  ErrorSignals  is  KxN  in  size  and  has  the  difference  between  Estimated 
and 

%  ToBeEstimated  signals 

%  MeanErrors  is  a  K  size  vector  that  has  the  mean  square  error  for  each 
to  be 

%  estimated  signal  in  it,  normalized  by  the  instantaneous  amplitude  of 
the  signal 

function  [EstTempCoef f ,  FullTempCoef f ,  Error]  = 

SensorBasedTempCoef f (XLoc,  YLoc,  XAxis,  YAxis, 

Pjj>DModes,  PODMean,  Snapshots) 

%  *[EstTempCoeff ,  FullTempCoef f ,  Error]  ... 

%  =  SensorBasedTempCoef f (XLoc,  YLoc,  XAxis,  YAxis, 

PODModes,  PODMean,  Snapshots) 

% 

%  Calculates  the  Temporal  Coefficients  based  on  both  the  full  flow  field 
%  contained  in  the  Snapshots, 

%  as  well  as  the  Flow  Field  sampled  at  the  XLoc,  YLoc  Locations 
%  Error  is  the  difference  between  the  two 

%  Uses  least  square  fitting  to  estimate  the  Temp  Coeff  for  both  full  flow 
%  field  and  sensor  based  estimation 

function  [ModeCoeff,  Indices]  =  POD_Model (TimeCoef f , 

MaxMode) 

%  Creates  a  low  order  model  dAn/dt  =  f(C0  +cliAi  +  c2ijAij  +c3ijkAijk) 
of  cubic  order 

%  Indices  returns  a  2D  Array,  the  columns  are  first  second  third 
coefficient  and  so  on 
%  the  rows  are  one  term  each 

%  Calculates  all  linear  and  quadratic  coefficients 
%  Modified  to  include  only  cubic  terms  up  to  Mode  Number  MaxMode 

function  [FilteredTime,  FilteredTempCoef f ]  = 

SplineFilterTempCoef f (Time,  TempCoeffs,  NPoints) 

% 

%  [FilteredTime,  FilteredTempCoef f]  =  SplineFilterTempCoef f (Timh , 
TempCoeffs,  NPoints) 

%  uses  spline_f ilter  to  filter  each  Temporal  Coefficent 

%  NPoints  determines  the  width  of  filter  kernel  in  data  points,  and  has 
one 

%  element  for  each  TempCoeff  to  allovf  different  filter  settings  for  each 
%  mode 

%  the  resulting  FilteredTempCoef f  are  shorter  in  time  by  two  times  the 
%  maximum  filter  kernel  width 

function  CFD_Controller_Main_V3 (CaseName, varargin) 

%  CFD_Controller_Main_V3 (CaseName)  j?Vs  before 

%  CFD_Controller_Main_V3 (CaseName, NodeNumber)  NodeNumber:  hostnode 

of  cobalt  run 
% 

%  interfaces  to  Cobalt  CFD  code  for  closed  loop  flow  control  ' 

%  filepath  is  the  fully  qualified  path  to  the  hdf  file  created  by  cobalt 
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%  Maxlteration  is  the  number  of  time  steps  for  this  run 

%  V3  supports  unstructured  POD  modes  and  feedback  controlled  Boundary 

% 

%  varargin  added  for  support  of  remote  execution 

function  [FileStructure]  =  Read_HDF_Cobalt (Filename, 

varargin) 

% 

%  [FileStructure]  =  Read__HDF_Cobalt  (Filename) 

%  [FileStructure]  =  Read_HDF_Cobalt (Filename,  Skip) 

%  [FileStructure]  =  Read_HDF_Cobalt (Filename ,  NStart ,NTimeSteps) 

%  [FileStructure]  =  Read_HDF_Cobalt (Filename,  Skip,  NStart ,NTimeSteps) 

%  Reads  all  the  non-empty  SDS'es  in  an  hdf  mtap  file  created  by  Cobalt  V3 
%  If  NStart  and  NTimeSteps  are  specified,  only  the  N  Time  steps  starting 
at  NStart 
%afiH  be  read 

%  -If  Skip  is  specified,  only  every  Skip-th  data  point  in  space  will  be 
%  loaded 

%  FileStructure:  Name  of  a  structure  that  will  hold  all  the  variables  in 
the  hdf  file 
%  c  2004  Stefan  Siegel 


t 
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